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1.  Research  Abstract 


A  unified  framework  is  proposed  for  providing  a  systematic  scheme  for  generating  the 
data  association  hypotheses  efficiently  in  the  target-oriented,  measurement-oriented,  and  track- 
oriented  approaches  to  multitarget  tracking.  A  fast  recursive  algorithm  for  computing  the  a 
posteriori  probabilities,  suitable  for  implementation  in  a  distributed  multiprocessor  system  is 
developed  and  its  links  to  the  theory  of  permanents  is  estabhshed.  An  analysis  of  this 
algorithm  reveals  its  superiority  over  existing  ones  in  the  average  case.  In  the  related 
problem  of  direction-of-arrival  estimation,  a  new  non-search-type  subspace  method,  called  the 
PESS  method,  is  proposed.  This  method  exploits  the  structure  of  the  steering  matrix  more 
thoroughly  to  yield  a  residual-error  theoretically  shown  to  be  either  less  than  or  equal  to  that 
obtained  by  LS-ESPRIT.  Furthermore,  simulation  conducted  on  several  sets  of  data  showed 
that  the  PESS  method  outperforms  the  TLS-ESPRIT  method.  Constraints  for  forcing  all  roots 
of  a  polynomial  to  the  unit  circle  are  obtained  for  more  reliable  estimation  especially  in  the 
low  SNR  case.  Finally,  for  improved  preprocessing  to  facilitate  tracking,  a  theoretical 
analysis  is  proposed  to  evaluate  the  robustness  of  a  TLS  algorithm,  developed  earlier,  for 
image  reconstruction  from  a  sequence  of  undersampled  noisy  and  blurred  frames. 
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2.  Research  Description 


(a)  Background 

To  solve  the  problem  of  data  association  between  targets  and  measurements,  three 
typical  approaches  have  been  reported  in  the  literature.  The  first  is  called  a  target-oriented 
approach  in  which  each  measurement  is  assumed  to  have  originated  from  either  a  known 
target  or  clutter,  as  in  the  Joint  Probabilistic  Data  Association  Filter  (JPDAF)  [1].  The  scond 
is  called  a  measurement-oriented  approach  in  which  each  measurement  is  hypothesized  to 
have  originated  from  either  a  known  target,  a  new  target,  or  clutter  [2].  The  third  approach  is 
referred  to  as  track-oriented  where  each  track  is  hypothesized  to  be  either  undetected, 
terminated,  associated  with  a  measurement,  or  linked  to  the  start  of  a  maneuver  [3]. 

In  recent  years,  a  lot  of  attention  has  been  given  to  the  task  of  improving  the 
computational  efficiency  of  a  multitarget  tracking  algorithm.  Fitzgerald  [4]  proposed  an  ad 

hoc  formula  to  approximately  compute  the  p‘’s  in  the  JPDAF.  For  j  >  Q,  pj  is  the  a 

posteriori  probability  that  measurement  j  originated  from  target  t  ■  Fitzgerald’s  formula 
breaks  down  when  four  targets  are  required  to  be  tracked  [5].  Sengupta  and  litis  [5] 
developed  an  analog  neural  network  to  emulate  the  JPDAF.  They  showed  that  the  neural 
network  is  capable  of  handling  two  to  six  targets  and  three  to  twenty  measurements. 

However,  the  heuristic  nature  of  their  approach  makes  implementation  difficult  [6]. 
Alternatively,  Nagarajan,  et  al,  [7]  arranged  the  hypotheses  in  the  measurement-oriented 
approach  [2]  in  a  special  order  so  that  the  probabilities  of  the  hypotheses  are  proportional  to 
the  product  of  certain  probability  factors  already  evaluated.  The  algorithm  locates  the  N 
globally  best  hypotheses  without  evaluating  all  of  them.  In  [8],  Fisher  and  Casasent 
developed  a  fast  joint  probabilistic  data  association  (JPDA)  algorithm.  In  their  algorithm,  the 

computation  of  the  p‘’s  was  implemented  using  enormous  number  of  vector  inner  product 

operations.  Since  the  vector  inner  product  operation  could  be  easily  realized  on  an  optical 
processor,  they  proposed  a  specialized  optical  processor  to  approximately  implement  the  fast 
JPDA  algorithm.  Recently,  Zhou  and  Bose  [9]  proposed  a  depth-first  search  (DFS)  approach 

to  efficiently  compute  the  pj’s  in  the  JPDAF.  Their  algorithm  requires  much  less 

computation  than  the  fast  JPDA  algorithm  in  the  average  case,  when  there  are,  on  average, 
two  to  three  measurements  inside  the  validation  gate  of  a  target.  Although  this  DFS 
algorithm  is  much  more  efficient  than  the  fast  JPDA  algorithm  in  the  average  case,  it  is 
specifically  directed  towards  implementation  through  a  system  with  a  powerful  centralized 
processor,  normally  used  in  ground-based  tracking  systems. 

(b)  Research  Contributions 

In  this  research,  an  efficient  algorithm  has  been  developed  to  compute  the  a  posteriori 
probabilities  of  the  origins  of  measurements  in  the  joint  probabilistic  data  association  filter 
(JPDAF).  The  inherited  parallelism  of  this  algorithm  enables  it  to  be  suitable  for 
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implementation  in  a  multiprocessor  system.  In  this  algorithm,  the  a  posteriori  probability  of 
the  origin  of  each  measurement  in  the  JPDAF  is  decomposed  into  two  parts.  The 
computation  of  one  part  becomes  trivial  and  the  algorithm  developed  here  is  implemented  on 
the  other  part,  which  is  shown  to  be  related  to  permanents.  The  computational  complexity  of 
this  algorithm  has  been  analyzed  in  the  worst  case  as  well  as  in  the  average  case.  It  has  been 
concluded  that  this  algorithm  is  more  efficient  than  other  existing  ones  in  the  average  case. 
The  results  are  fully  documented  in  the  very  recent  peer-reviewed  publication  [10]. 

Another  important  contribution  emerging  from  this  research  is  the  development  of  a 
unified  framework  to  provide  a  comprehensive  understanding  of  the  problem  of  data 
association  in  multitarget  tracking.  Specifically,  the  DFS  algorithm  in  [9],  which  was 
developed  under  the  sponsorship  of  the  grant  preceding  this  research,  has  been  adapted  for 
efficiently  generating  the  data  association  hypotheses  in  the  measurement-oriented  and  track- 
oriented  approaches,  where  the  total  number  of  data  association  hypotheses  is  expected  to 
increase  drastically  over  that  in  the  target-oriented  approach.  However,  reduction  in  the 
overall  computational  cost  may  be  feasible  from  observations  on  the  conditional  likelihood  of 
data  association  hypotheses.  In  the  target-oriented  approach,  the  conditional  likelihood  of 
each  data  association  hypothesis  is  unique.  When  targets  are  grouped  into  clusters,  this 
uniqueness  property  does  not  hold  for  the  measurement-oriented  and  track-oriented 
approaches.  Two  specialized  DFS  algorithms  which  can  efficiently  identify  the  data 
association  hypotheses  with  identical  conditional  likelihood  in  the  measurement-oriented  and 
track-oriented  approaches  have  been  developed  [11]. 

The  problems  of  direction-of-arrival  (DOA)  estimation  and  parameter  estimation  of 
sinusoids  in  noise  have  been  tackled  extensively  by  subspace-based  methods  during  the  last 
decade  ever  since  the  harmonic  retrieval  method  was  proposed  by  Pisarenko  in  1973  followed 
by  Schmidt’s  doctoral  dissertation  in  1979  which  led  to  a  formal  presentation  of  subspace- 
based  methods  in  the  open-literature  in  1985.  A  new  approach  to  parameter  estimation  based 
on  signal-parameter-selectivity  of  the  signal  subspace  (PESS)  was  proposed  recently  [12],[13]. 
This  method  exploits  the  strueture  of  the  steering  matrix  more  thoroughly.  The  residual  error 
is  theoretically  analyzed  and  it  has  been  shown  that  in  the  PESS  method  this  residual  error  is 
less  than  or  equal  to  that  of  the  LS-ESPRIT.  In  simulation,  it  is  shown  that  the  PESS  method 
outperforms  the  TLS-ESPRIT  method  [14]. 

In  DOA  estimation  and  the  related  problem  of  parameter  estimation  of  sinusoids  in 
noise,  the  need  for  constraining  the  coefficients  of  a  polynomial  so  that  its  roots  fall  on  the 
unit  circle  occurs.  While  simple  symmetry  conditions  suffice  in  a  lot  of  cases  [15],  the  need 
for  more  powerful  constraints  arises  in  low  SNR  situations.  Investigations  into  such 
constraints  lead  to  a  set  of  necessary  conditions  on  the  coefficients  of  a  polynomial  for  all  its 
roots  to  lie  on  a  unit  circle  [16].  The  powerful  mathematical  structure  built  around  the  theory 
of  resultants  used  in  multidimensional  systems  theory  for  a  variety  of  purposes  [17],  [18]  are 
again  used  to  attain  the  desired  objective  here. 

Finally,  to  link  image  processing  to  target  tracking  for  better  tracker  performance,  a 
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recursive  procedure  based  on  total  least  squares  (TLS)  theory  to  reconstruct  a  high  resolution 
image  from  a  sequence  of  low  resolution  noisy  and  blurred  frames  (incorporating  the  inherent 
uncertainty  in  displacement  estimations  of  the  frames  with  respect  to  a  reference  frame)  was 
developed  in  [19]  as  part  of  research  conducted  under  the  sponsorship  of  the  previous  ONR 
Grant  N00014-86-K-0542.  During  the  final  phase  of  this  research,  a  theoretical  analysis  was 
conducted  to  evaluate  the  robustness  of  the  TLS  algorithm  developed  in  [19].  It  was  shown 
that  with  certain  assumptions  on  the  noise,  the  image  reconstructed  using  the  TLS  algorithm 
has  minimum  variance  with  respect  to  all  unbiased  estimates.  Furthermore,  the  quality  of  the 
reconstructed  image  improved  with  increase  in  the  number  of  undersampled  frames.  In  the 
case  of  blurred  frames,  higher  resolution  images  may  be  reconstructed  using  the  TLS 
algorithm  with  post-deblurring  [20]. 
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163  Sensor  Array  Processing _ _ _ _ _ 

N.  K,  Bose  and  I.  H.  Sibul 

I  Multidimensional  signal  processing  tools  apply  to  aperture  and  sensor  array  processing.  Planar 
sensor  arrays  can  be  considered  to  be  sampled  apertures.  Three-dimensional  or  volumetric  arrays 
can  be  viewed  as  multidimensional  spatial  filters.  Therefore,  the  topics  of  sensor  array  processing, 
aperture  processing,  and  multidimensional  signal  processing  can  be  studied  under  a  unified  for¬ 
mat  The  basic  function  of  the  receiving  array  is  transduction  of  propagating  waves  in  the  medium 
into  electrical  signals.  Propagating  waves  are  fundamental  in  radar,  commumcation,  optics,  sonar, 
and  geophysics.  In  electromagnetic  applications,  basic  transducers  are  antennas  and  arrays  of 
antennas.  A  large  body  of  literature  that  exists  on  antennas  and  antenna  arrays  can  be  exploited  in 
the  areas  of  aperture  and  sensor  array  processing.  Much  of  the  antenna  literature  deals  with  trans¬ 
mitting  antennas  and  their  radiation  patterns.  Because  of  the  reciprocity  of  transmitting  and 
receiving  transducers,  key  results  that  have  been  developed  for  transmitters  can  be  used  for  analysis 
of  receiver  aperture  and/or  array  processing.  Transmitting  transducers  radiate  energy  in  desired 
directions,  whereas  receiving  apertures/arrays  act  as  spatial  filters  that  emphasize  signals  from  a 
desired  look  direction  while  discriminating  against  interferences  from  other  directions.  The  spatial 
filter  wavenumber  response  is  called  the  receiver  beam  pattern.  Transmitting  apertures  are  charac- 
^  terized  by  their  radiation  patterns. 

j  Conventional  beamforming  deals  with  the  design  of  fixed  beam  patterns  for  given  specifications. 

^  Optimum  beamforming  is  the  design  of  beam  patterns  to  meet  a  specified  optimization  criterion. 
It  can  be  compared  to  optimum  filtering,  detection,  and  estimation.  Adaptive  beamfonners  sense 
their  operating  environment  (for  example,  noise  covariance  matrix)  and  adjust  beamformer 
‘  parameters  so  that  their  performance  is  optimized  [Monzingo  and  Miller,  1980].  Adaptive  beam- 
formers  can  be  compared  with  adaptive  filters.  " 

Multidimensional  signal  processing  techniques  have  found  wide  application  in  seismology 
'vhere  a  group  of  identical  seismometers,  called  seismic  arrays,  are  used  for  event  location,  studies 
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of  the  earth’s  sedimentation  structure,  and  separation  of  coherent  signals  from  noise,  which  some¬ 
times  may  also  propagate  coherently  across  the  array  but  with  different  horizontal  velocities _ by 

employing  velocity  filtering  [Claerbout,  1976).  Velocity  filtering  is  performed  by  multidimen¬ 
sional  filters  and  allows  also  for  the  enhancement  of  signals  which  may  occupy  the  same  wavenum¬ 
ber  range  as  noise  or  undesired  signals  do.  In  a  broader  context,  beamforming  can  be  used  to  sepa- 
rate  signals  received  by  sensor  arrays  based  on  frequency,  wavenumber,  and  velocity  (speed  as  weQ 
as  direction)  of  propagation.  Both  the  transfer  and  unit  impulse-response  functions  of  a  velocity 
filter  are  two-dimensional  functions  in  the  case  of  one-dimensional  arrays.  The  transfer  function 
involves  frequency  and  wavenumber  (due  to  spatial  sampling  by  equally  spaced  sensors)  as  inde¬ 
pendent  variables,  whereas  the  unit  impulse  response  depends  upon  time  and  location  within  the 
array.  Two-dimensional  filtering  is  not  limited  to  velocity  filtering  by  means  of  seismic  array.  Two- 
dimensional  spatial  filters  are  frequently  used,  for  example,  in  the  interpretation  of  gravity  and 
magnetic  maps  to  differentiate  between  regional  and  local  features.  Input  data  for  these  filters 
may  be  observations  in  the  survey  of  an  area  conducted  over  a  planar  grid  over  the  earth  s  surface. 
Two-dimensional  wavenumber  digital  filtering  principles  are  useful  for  this  purpose.  Velocity  fil¬ 
tering  by  means  of  two-dimensional  arrays  may  be  accomplished  by  properly  shaping  a  three- 
dimensional  response  function  H(A:i,/:2,0)).  Velocity  filtering  by  three-dimensional  arrays  may  be 
accomplished  through  a  four-dimensional  function  H(ki,k2yk^yG>)  as  explained  in  the  following 
subsection. 


spatial  Arrays,  Beamformers,  and  FIR  Filters 

A  propagating  plane  wave,  s(x,f),  is,  in  general,  a  function  of  the  three-dimensional  space  variables 
and  the  time  variable  (xi^XiyX^)  4  x  and  the  time  variable  t.  The  4-D  Fourier  transform  of  the  sta¬ 
tionary  signal  s(x,r)  is 


..  S(k,(o)  -  J  j  J  J  (16.3) 

which  is  referred  to  as  the  wavenumber-frequency  spectrum  of  s(x,r),  and  (fci,^2>^3)  =  ^  denotes  the 
wavenumber  variables  in  radians  per  unit  distance  and  ©  is  the  frequency  variable  in  radians  per 
second.  If  c  denotes  the  velocity  of  propagation  of  the  plane  wave,  the  following  constraint  must  be 
satisfied  ..  .  >  : 


kl  +  kl  +  k\ 

>  '  ■  /.  ■/  if! 


If  the  4-D  Fourier  transform  of  the  unit  impulse  response  h{x,t)  of  a  4-D  linear  shift- invariant 
(LSI)  filter  is  denoted  by  H(  fc,©),  then  the  response  y(x,r)  of  the  filter  to  s(x,t)  is  the  4-D  linear  con* 
volution  of  hix,t)  and  5(x,r),  which  is,  uniquely,  characterized  by  its  4-D  Fourier  transform 


The  inverse  4-D  Fourier  transform,  which  forms  a  4-D  Fourier  transform  pair  with  Eq.  (16.3),  ^ 

1 


1 .70:.;. 


'  s(x.  t)  =  ^  r  f  r  f 


(16.3) 
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It  is  noted  that  S(k,0))  in  Eq.  (16.3)  is  product  separable,  i.e.,  expressible  in  the  form 
.  .  S(k,®)  =  S,(fci)S2(k2)S3(fc3)S4(CO) 
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linear  array 


P  FIGURE  16.23  Uniformly  weighted  linear  array. 

it  where  each  function  on  the  right-hand  side  is  a  univariate  function  of  the  respective  independent 
fe'Wiable,  if  and  only  if  5(x,f)  in  Eq.  (16.3)  is  also  product  separable.  In  beamfbrming,  SM  in 
%  Eq.  (16.6)  would  be  the  far-field  beam  pattern  of  a  linear  array  along  the  x.-axis.  For  example,  the 
i  normalized  beam  pattern  of  a  uniformly  weighted  (shaded)  linear  array  of  length  L  is 


kL  sin  6 


Sim  = 


where  X  =  {iKlk)  is  the  wavelength  of  the  propagating  plane  wave  and  0  is  the  angle  of  arrival  at 
array  site  as  shown  in  Fig.  16.23.  Note  that  9  is  explicitly  admitted  as  a  variable  in  S(fc,0)  to  allow  for 
the  possibility  that  for  a  fixed  wavenumber,  the  beam  pattern  could  be  plotted  as  a  function  of  the 
angle  of  arrival.  In  that  case,  when  0  is  zero,  the  wave  impinges  the  array  broadside  and  the  nor¬ 
malized  beam  pattern  evaluates  to  unity. 

The  counterpart,  in  aperture  and  sensor  array  processing,  of  the  use  of  window  functions  in 
spectral  analysis  for  reduction  of  sidelobes  is  the  use  of  aperture  shading.  In  aperture  shading,  one 
simply  multiplies  a  uniformly  weighted  aperture  by  the  shading  function.  The  resulting  beam  pat¬ 
tern  is,  then,  simply  the  convolution  of  the  beam  pattern  of  the  uniformly  shaded  volumetric  array 
and  the  beam  pattern  of  the  shading  function.  Fourier  transform  relationship  between  the  station¬ 
ary  signal  5(x,f)  and  the  wavenumber  frequency  spectrum  S(k,(D)  aUows  one  to  exploit  high-reso- 
lution  spectral  analysis  techniques  for  the  high-resolution  estimation  of  the  direction  of  arrival 
[Pillai,  1989].  The  superscript  t,  and  H  denote,  respectively,  complex  conjugate,  transpose,  and 
conjugate  transpose. 


Discrete  Arrays  for  Beamforming 

;  An  array  of  sensors  could  be  distributed  at  distinct  points  in  space  in  various  ways.  Line  arrays,  pla- 
I  nar  arrays,  and  volumetric  arrays  could  be  either  uniformly  spaced  or  nonuniformly  spaced, 
including  the  possibility  of  placing  sensors  randomly  according  to  some  probability  distribution 
^nnetion.  Uniform  spacing  along  each  coordinate  axis  permits  one  to  exploit  the  well-developed 
multidimensional  signal  processing  techniques  concerned  with  filter  design,  DFT  computation  via 
<  and  high-resolution  spectral  analysis  of  sampled  signals  [Dudgeon,  1977].  Nonuniform  spac- 
sometimes  might  be  useful  for  reducing  the  number  of  sensors,  which  otherwise  might  be  con- 
to  satisfy  a  maximum  spacing  between  uniformly  placed  sensors  to  avoid  grating  lobes 
^  ^  to  aliasing,  as  explained  later.  A  discrete  array,  uniforinly  spaced,  is  convenient  for  the  synthe- 
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sis  of  a  digital  filter  or  beamformer  by  the  performing  of  digital  signal  processing  operation^  1 1 
(namely  delay,  sum,  and  multiplication  or  weighting)  on  the  signal  received  by  a  collection  of  sen,  i  f 
sors  distributed  in  space.  The  sequence  of  the  nature  of  operations  dictates  the  types  of  beam 
former.  Common  beamforming  systems  are  of  the  straight  summation,  delay-and-sum,  and 
weighted  deiay-and-sum  types.  The  geometrical  distribution  of  sensors  and  the  weights  w,  assocj. 
ated  with  each  sensor  are  crucial  factors  in  the  shaping  of  the  filter  characteristics.  In  the  case  of  ^ 
linear  array  of  N  equispaced  sensors,  which  are  spaced  D  units  apart,  starting  at  the  origin 
Xi  =  0,  the  function 

N 


( 16,8) 


n=0 


becomes  the  array  pattern,  which  may  be  viewed  as  the  frequency  response  function  for  a  finite 
impulse  response  (FIR)  filter,  characterized  by  the  unit  impulse  response  sequence  { w„\.  In  the  case 
when  w„  =  1,  Eq.  (16.8)  simplifies  to 
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Sin 
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If  the  N  sensors  are  symmetrically  placed  on  both  sides  of  the  origin,  including  one  at  the  origin, 
and  the  sensor  weights  are  =  1,  then  the  linear  array  pattern  becomes 


Sin 


'’kiND^ 


W(k,)  = 


N 


sm 


For  planar  arrays,  dirert  generalizations  of  the  preceding  linear  array  results  can  be  obtained.  To 
wit,  if  the  sensors  with  unity  weights  are  located  at  coordinates  {kD,  ID),  where  ^  =  0,  ±  1,  ±2, . .  •» 
±[(N - 1)/2],  and  /  =  0,  ±  1,  ±2, ...» ±[(M- 1)/2],  for  odd  integer  values  of  N and  M,  then  the  arrar 
pattern  function  becomes 

W{ki,k2)  = -^  ^  '^eipi-jikikD  +  k^lD} 
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;  Routine  generaliz^ons  to  3-D  spatial  arrays  are  also  possible.  The  array  pattern  functions  for 
;  o4er  geometrical  istnbutions  may  also  be  routinely  generated.  For  example,  if  unit  weight  sen- 
,  jors  are  located  at  the  six  vertices  and  the  center  of  a  regular  hexagon,  each  of  whose  sides  is  D  units 

long,  then  the  array  pattern  fonction  can  be  shown  to  be 


1  +  2co,^,D+4cos^£cos^^ 
7  '  2  ■■  2- 


(16.11) 


J:  Ih,  y  patten,  toct.oa  .e.ola  how  selective  a  particja,  bca,nfonninga,„„,.i,,to*ecaae 

5  rf.  W.C.1  .nay  foncPoo  ahown  m  B,.  (16.9).  the  beamwidd.,  ridch  la  the  w  Ju.  of  dw  oiain  lobe 

;  •fthe.n.ypanerniainvenelypr^rtoaoJtod.eanayap.nore.Becaoaeoftheperiodidtyofth! 

. .  «iay  pattern  Ijinction.  the  mam  lobe  is  repeated  at  intervals  of  2ir/D.  Hiese  tepeddve  l^s  ate 
j  aDrf  grating  lobes,  whose  existence  may  be  interpreted  in  terms  of  spatial  fi^uency  aliasing 

;  .  ,<»6.irig  from  a  sampling  mrerval  D  due  to  die  N  receiving  sensors  located  a.  diLm^^^^^^ 

;  ^ce.Ifthe  spacing  D  between  sensors  satisfies  m 


(i6.12) 


I  .here  X  is  foe  sinallest  wavelength  component  in  the  signal  received  by  the  array  of  sensors  then 
^^ating  fobes  have  no  effect  on  the  received  signal.  A  plane  wave  of  unit  ampUtude  which  is  inci- 

■^■^.P-lucesou^udatthesen. 

I  s(0)4s9  =  [exp(;0)  exp(j^,Dsin  6) . . .  exp(jki(N-  l)Dsin  6)]'  (16.13) 

wavenumben  In  array  processing,  the  array  output  may  be  viewed  as  the 

I  ■W«mc.Iongadi,ecd„ncharac.c,i^byth.anglced.t,J„;:niS^ 


Be  -  (yrtdl.Sa)  =  ^  w.  exp(;l:,to  sin  6) 


(16.14) 


14  dons  lAhmed  anlferrosy  l7  M  f  despite  certiun  pemirba- 

K.-  h.  5  +  fh  I  f  component  s**  of  s*  to  belong  to  an  interval 

I  ».  w  wh  kh  Wilf  ““  *0  “dsdmce  of  at  least  one  weight  vec- 

^  '*ich  maybe  tackleH  h  *  •  Problem  can  be  translated  into  an  optimization  problem, 

I  be  tackled  by  minimizing  the  value  ofthe  array  output  power 


P(0)  =  w"(0)i?w(0) 


(16.15) 


iH«»l«tk„s’cdM  .fdve  "oise-corrupKd  signal  an.ocorreladon  matrix. 
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IS  called  the  minimum  variance  beamformer  and  is  given  by 
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^^^■responding  power  output  is 
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PuviQ)  = 


s"(0)i?*'s(9) 


The  minimum  variance  power  as  a  function  of  0  can  be  used  as  a  form  of  the  data- adaptive  esti 
mate  of  the  directional  power  spectrum.  However,  in  this  mode  of  solution,  the  coefficient  vector  i* 
unconstrained  except  at  the  steering  direction.  Consequently,  a  signal  tends  to  be  regarded  as  au 
unwanted  interference  and  is,  therefore,  suppressed  in  the  beamformed  output  unless  it  is  almost 
exactly  aligned  with  the  steering  direction.  Therefore,  it  is  desirable  to  broaden  the  signal  accep. 
tance  angle  while  at  the  same  time  preserving  the  optimum  beamformer  s  ability  to  reject  noise  and 
interference  outside  this  region  of  angles.  One  way  of  achieving  this  is  by  the  application  of  the 
principle  of  superdirectivity. 

Discrete  Arrays  and  Polynomials 

It  is  common  practice  to  relate  discrete  arrays  to  polynomials  for  array  synthesis  purposes  [Stein¬ 
berg,  1976].  For  volumetric  equispaced  arrays  (it  is  only  necessary  that  the  spacing  be  uniform 
along  each  coordinate  axis  so  that  the  spatial  sampling  periods  D/  and  Dj  along,  respectively,  the  rth 
and  ;th  coordinate  axes  could  be  different  for  i  ^  ;“),  the  weight  associated  with  sensors  located  at 
coordinate  is  denoted  by  wihyhyh)’  The  function  in  the  complex  variables 

and  Z3)  that  is  associated  with  the  sequence  {w(*i*J2>*3)}  is  the  generating  function  for  the  sequence 
and  is  denoted  by 


In  the  electrical  engineering  and  geophysics  literature,  the  generating  function  Mzj, 73,73)  is  some¬ 
times  called  the  z-transform  of  the  sequence  When  there  are  a  finite  number  of  sensors, 

a  realistic  assumption  for  any  physical  discrete  array,  W( 71,2^,73)  becomes  a  trivariate  polynomial 
In  the  special  case  when  Mh^hyb)  is  product  separable,  the  polynomial  Wizi^ZiyZ^)  is  also  product 
separable.  Particularly,  this  separability  property  holds  when  the  shading  is  uniform,  i.e., 

-  1.  When  the  support  of  the  uniform  shading  function  is  defined  by  ii  =  0,l,...,Ni-l,i:  = 
0, 1, . . . ,  N2-  1,  and  I3  =  0,1, . . . ,  N3-  1,  the  associated  polynomial  becomes 

N,-l  N,-i  N3-.1  3  , 

W{Zi,  Z2,  Zj)  =  =  H  ■  ■  _  (16.19) 

Ii  =  0  12  =  0  iy  =  0  i=l 

In  this  case,  all  results  developed  for  the  synthesis  of  linear  arrays  become  directly  applicable  to  the 
synthesis  of  volumetric  arrays.  For  a  linear  uniform  discrete  array  composed  of  N  sensors  with 
intersensor  spacing  Di  starting  at  the  origin  and  receiving  a  signal  at  a  known  fixed  wavenumber 
at  a  receiving  angle  0,  the  far-field  beam  pattern 

.  N-l 

S(it,,0)^S(e)= 

r=0 

may  be  associated  with  a  polynomial  2^o7i,  by  setting  Zi  =  This  polynomial  has  all  its 

zeros  on  the  unit  circle  in  the  7i -plane.  If  the  array  just  considered  is  not  uniform  but  has  a  weight¬ 
ing  factor  w^,  for  r  =  0, 1, . . . ,  Ni  -  1,  the  space  factor, 
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P’  again  be  associated  with  a  polynomial  Z'^LoVz^  By  the  pattern  ihultiplication  theorem,  it  is 
ssible  to  get  the  polynomial  associated  with  the  total  beam  pattern  of  an  array  with  weigh  ted  sen- 
^ts  by  multiplying  the  polynomials  associated  with  the  array  element  pattern  and  the  polynomial 
f  Related  with  the  space  factor  0(0).  The  array  factor  |Q(0)|^  may  also  be  associated  with  the  poly- 

'nomial  spectral  factor 


Ni-l  --w,-! 

\m\ 


(16.20) 


^ere  the  weighting  (shading)  factor,  is  allowed  to  be  complex.  Uniformly  distributed  apertures 
S  and  uniformly  spaced  volumetric  arrays  which  adimt  product  separable  sensor  weightings  can  be 
'Wed  by  using  the  well-developed  theory  of  linear  discrete  arrays  and  Aeir  associated  polyno- 
^  'ipial  When  the  product  separability  prppeirty  does  not  hold,  scopes  exist  for  applying  results  from 
^  multidimensional  systems  theory  [Bose,  1982]  concerning  multivariate  polynomials  to  the  synthe- 
j^ieproblem  of  volumetric  arrays.  ^ 


Velocity  FUtering  / 

Combination  of  individual  sensor  outputs  in  a  more  sophisticated  way  than  the  delay-and-sum 
technique  leads  to  the  design  of  multichannel  velocity  filters  for  linear  and  plan^  as  well  ^  spatial 
'  arrays.  Consider,  first,  a  linear  (1-D)  iriray  of  sensors,  which  will  be  used  to  implement  velocity  dis¬ 
crimination.  The  pass  and  rejection  zones  are  defined  by  straight  lines  in  the  (fci,(o) -plane,  where 

®  -  CO 

\  V  (v/sin  0) 

is  the  wavenumber,  (0  the  angular  frequency  in  radians/second,  V  the  apparent  velocity  on  the 
earth*s  surface  along  the  array  line,  v  the  velocity  of  wave  propagation,  and  0  the  horizontal  arrival 
direction.  The  transfer  function 


Hidiyk,)  = 


1. 

V  V 

0.  otherwise 


of  a  “pie-slice”  or  “fan”  velocity  filter  [Bose,  1985]  rejects  totally  wavenumbers  outside  the  range 
“iu/ V S  it|  <  |(o|/ Vand  passes  completely  wavenumbers  defined  within  that  range.  Thus,  the  trans¬ 
fer  function  defines  a  high-pass  filter  which  passes  signals  with  apparent  velocities  of  magnitude 
peater  than  V  at  a  fixed  frequency  (0.  If  the  equispaced  sensors  are  D  units  apart,  the  spatial  sam- 
results  in  a  periodic  wavemunber  response  with  period  =  1/ (2D).  Therefore,  for  a  specified 
*Pi"rent  velocity  V,  the  resolvable  wavenumber  and  frequency  bands  are,  respectively,  — 1/(2D)  ^  ki 
S  l/(2D)  and -V/(2D)  <  (O  <  V/(2D)  where  a)/(2D)  represents  the  folding  frequency  in  radians/ 
oecond. 

linear  arrays  are  subject  to  the  limitation  that  the  source  is  required  to  be  located  on  the 
^i.^®®ded  line  of  sensors  so  that  plane  wavefronts  approaching  the  array  site  at  a  particular  velocity 
the  individual  sensors,  assumed  equispaced,  at  arrival  times  which  are  also  equispaced.  In 
^fefemology,  the  equispaced  interval  between  successive  sensor  arrival  times  is  caUed  a  move-out  or 
•'P-out  and  equals  (D  sin  6)/v=  D/V.  However,  when  the  sensor-to-source  azimuth  varies,  two  or 
independent  signal  move-outs  may  be  present  Planar  (2-D)  arrays  are  then  required  to  dis- 
^®nnate  between  velocities  as  well  as  azimuth.  Spatial  (3-D)  arrays  provide  additional  scope  to 
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the  enhancement  of  discriminating  capabilities  when  sensor/source  locations  are  arbitrary.  In  such 
cases,  an  array  origin  is  chosen  and  the  mth  sensor  location  is  denoted  by  a  vector  (^i  m  -^3 '  and 
the  frequency  wavenumber  response  of  an  array  of  sensors  is  given  by 

1  N 

H(q),  =  —  y  HJ(o)  exp 

where  denotes  the  frequency  response  of  a  filter  associated  with  the  mth  recording  device 

(sensor).  The  sum  of  all  AT  filters  provides  flat  frequency  response  so  that  waveforms  arriving  from 
the  estimated  directions  of  arrival  at  estimated  velocities  are  passed  undistorted  and  other  wave¬ 
forms  are  suppressed.  In  the  planar  specialization,  the  2-D  array  of  sensors  leads  to  the  theorv  of 
3-D  filtering  involving  a  transfer  function  in  the  frequency  wavenumber  variables/  Jti,  and  k..  The 
basic  design  equations  for  the  optimum,  in  the  least-mean-square  error  sense,  frequeno’ 
wavenumber  filters  have  been  developed  [Burg,  1964].  This  procedure  of  Burg  can  be  routinelv 
generalized  to  the  4-D  filtering  problem  mentioned  above. 
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Defining  Terms 

Array  pattern:  Fourier  transform  of  the  receiver  weighting  function  taking  into  account  the  posi¬ 
tions  of  the  receivers. 

Beamformers:  Systems  commonly  used  for  detecting  and  isolating  signals  that  are  propagating 
in  a  particular  direction.  * 

Grating  lobes:  Repeated  main  lobes  in  the  array  pattern  interpretable  in  terms  of  spatial  fre¬ 
quency  aliasing. 

Velocity  filtering:  Means  for  discriminating  signals  from  noise  or  other  undesired  signals 
because  of  their  different  apparent  velocities.  f  ~  - 

Wavenumben  2Tr  (spatial  frequency  in  cycles  per  unit  distance). 


References 


K.M.  Ahmed  and  R.  J.  Evans,  “Robust  signal  and  array  processing,”  lEE  ProceedingSy  F:  Communi- 
cationsy  Radar,  and  Signal  Processing,  vol.  129,  no.  4,  pp.  297-302, 1982. 

N.  K.  Bose,  Applied  Multidimensional  Systems  Theory,  New  York:  Van  Nostrand  Reinhold,  1982. 
N.K.  Bose,  Digital  Filters,  New  York:  Elsevier  Science  North-HoUand,  1985.  Reprint  ed.,  Malabar, 
Fla.:  Krieger  Publishing,  1993. 

J.  P.  Burg,  “Three-dimensional  filtering  with  an  array  of  seismometers,”  Geophysics,  vol.  23,  no.  5, 
pp. 693-713, 1964.  s-  .  •  ..r 

J.  F.  Claerbout,  Fundamentals  of  Geophysical  Data  Processing,  New  York:  McGraw-Hill,  1976. 

D.  E.  Dudgeon,  “Fundamentals  of  digital  array  processing,”  Proc,  IEEE,  \oL  65,  pp.  898-904, 1972- 

R.  A.  Monzingo  and  T.  W.  Miller,  Introduction  to  Adaptive  Arrays,  New  York:  Wiley,  1 980. 

S.  M.  Pillai,  Array  Signal  Processing,  New  York:  Springer- Verlag,  1989. 

B.  D.  Steinberg,  Principles  of  Aperture  and  Array  System  Design,  New  York:  Wiley,  1976. 


An  Efficient  Algorithm 
Data  Association  in  Multitarget 
Tracking  V 

B.  ZHOU,  Member,  IEEE 

N.  K-  BOSE,  Fellow,  IEEE 
Tlie  Pennsylvania  State  University 


An  effldent  algorithm  is  developed  to  compute  the  a  posteriori 
probabiliUes  of  the  origins  of  measurements  in  the  joint 
probabilistic  data  association  Alter  (JPDAF).  The  inherited 
parallelism  of  this  algorithm  enables  it  to  be  suitable  for 
ingilementation  in  a  multiprocessor  system.  In  this  algorithm,  the 
a  posteriori  probability  of  the  origin  of  each  measurement  in  the 
JPDAF  is  decomposed  Into  two  parts.  The  computation  of  one  part 
becomes  trivial  and  the  algorithm  developed  here  is  implemented 
on  the  other  part,  which  is  shown  to  be  related  to  permanents. 

The  computational  complexity  of  this  algorithm  is  analyzed  in 
the  worst  case  as  well  as  in  the  average  case.  An  analysis  of  this 
algorithm  enables  us  to  conclude  that  this  algorithm  is  more 
efficient  than  other  existing  ones  in  the  average  case. 
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I.  INTRODUCTION 

In  a  cluttered  environment,  the  received 
measurements  may  not  all  arise  from  the  targets  of 
interest  Some  of  them  may  be  from  clutter  or  false 
alarm.  As  a  result,  there  always  exist  ambiguities  in 
the  association  between  the  previous  known  targets 
and  measurements.  To  solve  the  problem  of  data 
association  between  targets  and  measurements,  three 
typical  approaches  have  been  reported  in  the  literature. 
The  first  is  called  a  target-oriented  approach  in  which 
each  measurement  is  assumed  to  have  originated 
from  either  a  known  target  or  clutter,  as  in  the  joint 
probabilistic  data  association  filter  (JPDAF)  [1,  2]. 

The  second  is  called  a  measurement-oriented  approach 
in  which  each  measurement  is  hypothesized  to  have 
originated  from  either  a  known  target,  a  new  target, 
or  clutter  [3].  The  third  approach  is  referred  to  as 
track-oriented  where  each  track  is  hypothesized  to 
be  either  undetected,  terminated,  associated  with  a 
measurement,  or  linked  to  the  start  of  a  maneuver 
[4-6].  In  these  approaches,  the  number  of  data 
association  hypotheses  could  inaease  rapidly  with  the 
increase  in  the  number  of  targets  and  the  number  of 
measurements.  Therefore,  in  a  multitarget  tracking 
algorithm,  the  computational  cost  in  generating 
the  data  association  hypotheses  would  be  excessive 
when  the  number  of  targets  and  the  number  of 
measurements  are  large. 

In  recent  years,  a  lot  of  attention  has  been  given 
to  the  task  of  improving  the  computational  efficiency 
of  a  multitarget  tracking  algorithm.  Htzgerald  [7] 
proposed  an  ad  hoc  formula  to  approximately  compute 
the  /3's  in  the  JPDAF.  For  j  >  0,  is  the  a  posteriori 
probability  that  measurement  j  originated  from  target 
t.  For  j  =  0,  /3o  is  the  a  posteriori  probability  that  no 
measurement  originated  from  target  t.  Fitzgerald’s 
formula  breaks  down  when  four  targets  are  required 
to  be  tracked  [8].  Sengupta  and  litis  [8,  9]  developed 
an  analog  neural  network  to  emulate  the  JPDAF.  They 
showed  that  the  neural  network  is  capable  of  handling 
two  to  six  targets  and  three  to  twenty  measurements. 
However,  the  heuristic  nature  of  their  approach 
makes  implementation  difficult  [10].  Alternatively, 
Nagarajan,  et  al  [11],  arranged  the  hypotheses  in 
the  measurement-oriented  approach  [3]  in  a  special 
order  so  that  the  probabilities  of  the  hypotheses  are 
proportional  to  the  product  of  certain  probability 
factors  already  evaluated.  The  algorithm  locates  the 
N  globally  best  hypotheses  without  evaluating  all 
of  them.  In  [12],  Fisher  and  Casasent  developed 
a  fast  joint  probabilistic  data  association  (JPDA) 
algorithm.  In  their  algorithm,  the  computation  of 
the  was  implemented  using  enormous  number 
of  vector  inner  product  operations.  Since  the  vector 
inner  product  operation  could  be  easily  realized  on 
an  optical  processor,  they  proposed  a  specialized 
optical  processor  to  approximately  implement  the 
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fast  JPDA  algorithm.  Recently,  Zhou  and  Bose  [13] 
proposed  a  depth-first  search  (DFS)  approach  to  * 
efficiently  compute  the  ^s  in  the  JPDAF.  Their  ■* 
algorithm  requires  much  less  computation  than  the 
fast  JPDA  algorithm  in  the  average  case,  when  there 
are,  on  average,  two  to  three  measurements  inside ' 
the  validation  gate  of  a  target.  The  validation  gate  of 
a  target  is  a  region  with  a  prescribed  size  around  the 
predicted  position  of  the  target. 

As  shown  in  [13,  14],  the  performance  of  an 
approximation  of  the  JPDAF  degrades  drastically 
when  the  density  of  targets  is  high.  Therefore,  it  is 
important  to  implement  the  JPDAF  accurately  in 
a  dense  target  environment  Among  the  algorithms 
mentioned  above,  only  the  two  algorithms  proposed  in 
[12,  13]  have  the  potential  for  accurately  computing 
the  /3Js  in  the  JPDAF.  The  fast  JPDA  algorithm  in 
[12]  is  suitable  for  implementation  in  a  multiprocessor 
system.  Nevertheless,  the  computational  cost  is  very 
high  when  the  number  of  targets  and  the  number 
of  measurements  are  large  since  the  design  of  the 
algorithm  is  based  on  the  worst  case  scenario.  In 
the  worst  case,  all  measurements  fall  inside  the 
intersection  of  the  validation  gates  of  all  targets. 
Although  the  DFS  algorithm  [13]  is  much  more 
efficient  than  the  fast  JPDA  algorithm  in  the  average 
case,  it  is  specifically  directed  towards  implementation 
through  a  system  with  a  powerful  centralized 
processor,  normally  used  in  ground-based  tracking 
systems.  Our  objective  here  is  to  propose  an  algorithm 
which  not  only  performs  better  than  the  DFS  algorithm 
in  the  average  case,  but  which  is  also  suitable  for 
implementation  in  a  tracking  system  with  several 
spatially  distributed  microprocessors.  . 

A  brief  review  of  the  JPDAF  is  given  in  Section  II, 
followed  by  the  problem  formulation  in  Section  III. 

The  new  algorit^  is  developed  in  Section  IV.  In  - 
Section  V,  the  computational  complexity  of  the  new 
algorithm  is  analyzed  in  the  worst  case  as  well  as  in 
the  average  case.  Finally,  the  advantages  of  the  new 
algorithm  are  summarized  in  Section  VI. 

II.  REVIEW  OF  THE  jPDAF 

We  assume  that  there  are  n  targets  being  tracked 
at  time  index  k.  Let  the  dynamic  model  for  target  t  be 
described  by 


identically  distributed  Gaussian  prbces^^w^inowi'’' 
covariances:  ^  ^ 

£{v'(A:)(v'0'))'}  =  (4) 


where  6kj  =  I  ^  Ic  -  j  and  6k j  =  0  otherwise.  The  , 
prime  superscript  indicates  matrix  transpositibiL  The 
initial  target  state  x^(0)  is  assumed  to  be  riorm^y 
distributed  with  mean  J&^(0 1 0)  and  covariance 
P^(0  I  0).  It  is  also  assumed  to  be  independent  of  w^(k) 
and  for  all  A:  >  0. 

Suppose  that  m  measurements  are  received  at 
time  index  k.  In  a  cluttered  environment,  m  does 
not  necessarily  equal  n  and  it  may  be  difficult  to 
distinguish  whether  a  measurement  originated  from 
a  target  or  from  clutter.  It  is  reasonable  to  denote  a 
measurement  at  time  index  k  by 


if  z{k)  is  from  target  t 
if  z(/:)  is  from  clutter. 


(5) 


The  measurement  Tf{k)  is  usually  assumed  to  be 
uniformly  distributed  in  the  surveillance  region  and 
the  number  of  clutter  is  subject  either  to  Poisson 
distribution  or  to  uniform  distribution  [1].  In  this 
paper,  Poisson  distribution  is  assumed  for  the  purpose. 

In  the  JPDAF,  the  a  posteriori  probability  is 
computed.  For  j  >  0,  /3j  is  the  a  posteriori  probability 
that  validated  measurement  j  originated  fi-om  target  t. 
Pfj  is  the  a  posteriori  probab^ty  that  no  measurement 
originated  from  target  t.  A  validated  measurement  is 
one  which  is  either  inside  or  on  the  boundary  of  the 
vahdation  gate  of  a  target  Mathematically,  a  validation 
gate  is  defined  by 

(z(k)  -  2‘(k))'S‘{k)-\z(k)  -  2‘(k))  <  (6) 


where  z‘{k)  is  the  predicted  value  of  z(*)  for  target  t. 
The  error,  (z(k)- 2‘{k)),  is  the  iimovation  generated 
fi-om  z(k)  for  target  t,  S‘{k)  is  its  covariance  matrix, 
and  g  is  a  selected  threshold.  As  discussed  in  [2, 15], 
the  choice  g  >  VW  +  2  ensures  that  the  correct 
measurements  will  lie  within  the  gate  with  probabihty 
0.999.  The  dimension  of  z(k)  is  M‘.  The  inequality 
given  in  (6)  is  said  to  generate  a  validation  test  The 
result  of  the  test  is  kept  in  what  is  called  a  validation 
matrix  fi.  This  validation  matrix  f2  is  a  m  x  (n  -F 1) 
rectangular  matrix  defined  as  [2], 


x‘(k  1)  =  F‘{k)x\k)  +  G\k)w‘{k)  -1  (1) 

z‘{k)^H‘{k)x\k)  +  v\k),  t  = 

where  x‘{k)  is  the  iV' -dimensional  state  vector,  z*{k) 
is  the  M' -dimensional  measurement  vector  with  M' 
actually  independent  of  t,  F\k),  G‘(k')  and 
are  known  model  matrices,  w‘{k)  and  v‘(k)  are, 
respectively,  the  p^  and  Af^-dunensional  noise  vectors, 
which  are  assumed  to  be  zero-mean  independent 
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:  /=  1,2,..., m,  Uj,  =  1  if  .  ,5: 

measurement  j  is  inside  the  validation  gate  of  target  ^ 
cj}y=  for  /  =  and  r  = 

on  the  validation  matrix  fi,  data 
;  ^assoaaribn^hypotheses  (or  feasible  events  5s  [2])  are 
fgcmcrated  subject  to  the  following  two  restrictions: 

1)  Each  measurement  can  have  only  one  origin 
(cither  a  s^^cific  target  or  clutter). 

2)  No  more  than  one  measurement  originates  from 


a  target 

;  >s } 


This  leads  to  a  combinatorial  problem  where  the 
niimber  of  data  association  hypotheses  increases 
exponentially  with  n  and  m. 

'  Each  feasible  event  S  is  represented  by  a 
hypothesis  matrix  Cl  in  [2].  Cl  has  the  same  size  as  the 
validation  matrix  Q.  A  typical  element  in  is  denoted 
by  a)y7  where 


'1 

c2yf  =  I  1 


lo 


if  Wy,  =  1  and  measurement  j 
is  hypothesized  to  be  from 
clutter  for  t  =  0 

ifwy,  =  1  and  measurement  j  (8) 
is  hypothesized  to  be  associated 
with  target  r  for  f  0 

otherwise. 


In  (12),  H'C/c)  is  the  observation  matrix  for  target 
t  as  defined  in  (2)  and  fk  —  1)  is  the  one-step 
prediction  of  the  state  estimate  of  target  r.  It  is 
understood  that  the  normalizing  constant  c  in  (9) 
is  obtained  by  the  summation,  (ft)  P(S(p)  \  2). 
Therefore,  c  is  omitted  hereafter.  The  a  posteriori 
probability  /3‘j  is  computed  from  the  conditional 
probabilities  m  (9)  by 

£((1) 

«  =  I 

i=i  j  (1. 

for  j  =  l,2,...,m,  and  t  = 

The  equations  of  the  JPDAF  are  the  same  as  thos> 
of  a  standard  Kalman  filter  with  the  following  two 
changes. 

1)  The  innovation  vector  is  generated  by  using 

relation:  ^ 

z'(k)  =  ^P‘i(k)z‘j(k).  (1: 

;=i 

2)  The  update  equation  for  the  covariance  matrix 
P‘(k,k  I  k)  is  changed  to 


Corresponding  to  the  two  restrictions  above,  in  ll, 
there  is  at  most  one  unit  element  in  each  column 
except  for  r  =  0  and  there  is  exactly  one  unit  element 
in  each  row. 

After  each  Cl  is  obtained,  the  conditional 
probability  of  the  corresponding  data  association 
hypothesis  or  feasible  event  is  calculated  by  a  fonnula 
given  in  [2].  A  simplified  version  of  this  formula  is 
given  as 


P(£(^)  I Z)  =  P‘j  ■  .  ?•!':  V 

r.Oj.=i  ‘ 

for  j  =  l,2,...,m,  and  t  =  1,2,.'.’,^',“  ' 


where  Z  is  the  set  of  all  measurements  received  up  to 
current  time  index  c  is  a  normalizing  constant,  otj" 
is  the  number  of  targets  detected  in  this  feasible  ewht 
£,u)j,  =  1  indicates  that  measurement  j  is  associated  ' 
with  target  t  in  the  event,  and 


_  (N(z)-AS‘)Pd 

'  lo 


if  03  jt  -  1 
otherwise 


'(10) 


P^  =  A(l-Po)  =  Po 

for  y  =  l,2,...,m,  and  r  =  l,2,...,/i.  ^ 

In  (10)  and  (11),  A  is  the  clutter  density,  Po  is  the 
probability  of  detection,  and  N is  a  normal 
density  function  having  zero  mean  and  a  covariance 
matrix  with 


z](k)  =  Zj(k)  -  H‘  (k)Ji‘  (k  Ik-  1).  (12) 


P‘(k,k\k)  =  p*o(k)P‘(k,k\k-l) 

+  a-m))P‘*(k,k\k) 

+  K\k)W\k)K\k)'  (16 

where  K‘(k)  is  the  Kalman  gain  for  target  t, 

P‘*(k,k  I k)  =  [I-K‘(k)H‘(k)]P‘(k,k  \k-l) 

'  ,  (!"■ 

denotes  the  covariance  update  for  a  single  correct 
return,  and 

m 

-  z‘(k)z‘(ky.  (18) 
M 

In  general,  a  covariance  matrix  P‘(p,q\r)  is  defined 
as 

P‘(P>q  I  r)  =  E{[x‘(p)  -  i‘(p  i  r)][x\q)  -  3t‘{q  \  r)]'} 

(19) 

where  k‘(jplr)  =  E{x(j>)  |  measurements  received  up 
to  time  index  r}. 

In  the  next  seaion,  the  computation  of  the  /3jS  is 
reformulated  in  such  a  way  that  the  construction  of  the 
hypothesis  matrices  Cl  is  not  necessary. 

III.  PROBLEM  FORMULATION 

In  the  original  JPDAF  [2],  the  computation  of  /9  ■ 
consists  of  the  following  three  steps. 
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1)  Construct  a  validation  matrix  Q  and  compute  Pj 
using  (10)  and  (11). 

2)  Generate  data  association  hypotheses  and 
compute  P(£{Cl)  \  Z)  using  (9). 

3)  Compute  /3y  using  (13)  and  (14). 


In  this  section,  each  is  decomposed  into  two 
parts  so  that  its  computations  does  not  require  the 
generation  of  the  data  association  hypotheses  or  the 
construction  of  the  hypothesis  matrices  as  in  the 
original  JPDAF  [2]  and  the  DFS  algorithms  [13]. 

Let  A  denote  an  ordered  set  of  nonzero  integers 
associated  with  the  targets.  Let  B  denote  another 
ordered  set  of  non-zero  integers  which  are  the  indices 
of  measurements.  The  symbol  0  is  the  index  for  clutter 
Suppose  that  there  are  n  targets  and  m  measurements 
are  received  at  a  certain  time  instant  Then, 


B  =  {0,1,. .../n}. 

With  the  above  notation,  (13)  may  be  rewritten  as 


s(Ci) 


(22) 


conditioned  on  the  hypothesfei-  ti^ij^^ii.'tpw^i^'f  'y  is  Ic- 
associated  with  target  t. 


£s  are  mutuaUy  exclusive,  (22)  heieaifrW^  as 
beiov^  -  -•  ' 


£Cn)  '  ■  . 

=p  u  :  ■ 


=  p\£i,(Qj,)^  {  U  ^(ft)l|Z 


=  P(£j,(Wj,)\Z) 


f  \ 

\ 

(20) 

xP  j 

U 

nf/r(W/r)  \£j,(<^j,),Z 

(21) 

^£(n) 

J 

{ 


=  P{_£„{Co^,)  I  Z')P 


(J  £{Cl'^j,))\£i,{w„),Z 


for  j  =  l,2,...,m,  and  r  =  l,2,...,n. 

Since  Pj  is  a  factor  of  the  P(5(^l)  |  Z)s  when  ljji  =  1, 

I3‘j  may  be  decomposed  as 

l}'i(A.B)  =  P'jF(A\,.S\j)  ■  (23) 

for  /  =  l,2,...,m,  and  r  =  l,2,.,.,n 

,  ,  '  -■ 

where  the  difference  set  C\l  is  the  set  derived  from;iH>  ; 
C  by  deleting  element  /.  In  (23),  F(A\tyB\j)  is  a  fr^-  -  r 
function  of  P[  s  for  /  7^  j  and  r  ^t.  An  interpretation 
of  F(A\t,B\j)  may  be  obtained  using  the  definition  of 
P(P(i))  I  2)  [2], 

For  j  =  1,2,. ..,m,  and  t  =  0,1,2,. ..,n,  let  P),  ' 
denote  the  hypothesis  that  measurement  j  originated 
either  from  target  r  if  r  ^  0  or  from  clutter  if  t=  Q.,  _  : 
Then,  each  feasible  event  S  may  be  written  as  [2] 

^(^)='n^/>(0(^;>(0)  '  ^(24) 

^  ••  V 

where  p(t)  denotes  an  element  of  the  set  of 
permutations  of  {r}  ={0,1,2,.. .,n}.  Furthermore, 
let  be  the  hypothesis  matrix  obtained  from  Cl 
after  eliminating  the  yth  row  and  the  rth  column..  . 
P(£jt((^jt)  I  ^)  is  the  probability  of  the  hypothesis  ? 
that  measurement  j  originated  from  target  t  subject  ?  f 
to  the  assumption  that  data  association  between  the  - 
remaining  measurements  and  remaining  targets  is  I'- 
ignored,  and  I  (“;<). -2)  denotes 

the  probability  of  all  the  data  association  hypotheses 
among  the  remaining  measurements  and  targets  f  * 


=  P(£j,(u>i,)\Z)  Y, 

for  ;■  =  l,2,...,m,  and  t  =  l,2,...,n.  (25) 

Comparing  (25)  with  (23),  we  have 

P‘j=P(£j,(Qj,)\Z)  (26) 

F(A\t,B\j)=  P(^(%))|^/(W;v),2).  (27) 

Therefore,  F(A\t,B\j)  in  (23)  contains  the  sum  of 
the  conditional  probabilities  of  the  data  association 
hypotheses  among  all  the  targets,  the  measurements, 
and  clutter,  given  that  measurement  j  is  associated 
with  target  t.  In  the  case  of  the  probabilistic  data 
association  filter  (PDAF)  [15],  where  Pj  differs  from 
the  a  posteriori  probability,  /3j,  by  a  normalization 
constant,  F(A\tfB\j)  may  be  considered  to  be  due 
to  the  interference  from  other  targets  on  (23)  is 
compared  with  its  counterpart  in  the  PDAE 
Similarly,  /3J  may  be  decomposed  as 

P‘o  =  PoF(A\t,B),  for  r  =  l,2,...,n  (28) 

\  ^  ■^.5 

where  F(A\t,B)  contains  the  sum  of  the  conditional 
probabilities  of  the  data  association  hypotheses  among 
all  the  targets,  the  measurements,  and  clutter,  given 
that  target  t  is  not  associated  with  any  measurement 
In  (28),  F(A\t,B)  may  also  be  considered  as  the 
interference  on  /Sg  from  other  targets.  The  computation 
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^Street  a  validation  matrix  and  compute  Fj 
using  (lb)  and  (11). 

2)  Compute  F(^\r,5\y)  and  F(>1V,B). 

3)  Compute  /3j  using  (23)  and  (28). 

Tb  compute  for  t  =  using  (14)  instead 

of  (28),  the  p‘jS  need  to  be  normalized.  The  /3js 
computed  using  (23)  and  (28)  are  not  normalized.  The 
normalizing  constant  c  is  determined  by  the  following 
relationship. 

m 

'Y^j3‘j(A,B)  =  c,  for  r  =  l,2,...,/i.  (29) 

y=o 

Suppose  that  c  is  computed  in  the  t  =  1  case  as 

f;/3)(A5)  =  c.  (30) 

;=o 

The  same  constant  also  applies,  in  all  related 
cases  where  r  7^  1.  Then,  for  ;  =  0,l,...,m,  are 
normalized  as  below 

p)  - 1/3).  (31) 

Using  c,  for  t  =  2, and  }  =  can  also 

be  normalized  as  the  way  the  /3js  do.  Now,  Pq  for 
r  =  2,...,n,  are  ready  to  be  computed  using  (14).  In 
the  next  section,  an  algorithm  is  developed  to  compute 
F(A\t,B\j)  for  t  =  1,. and  ;  =  l,...,m,  and  ' 
F(A\1,B). 

IV.  ALGORITHM  DEVELOPMENT 

From  the  notation  introduced  in  the  last 
section,  F(A,B)  denotes  the  sum  of  the  conditional 
probabilities  of  the  data  association  hypotheses  among 
the  targets  in  A  and  the  measurements  in  B.  Using 
(22)  and  (27),  we  have  - 

•  f(A,B)  =  Y,P(^(^)\z)  '  ' 

. 


After  substituting  (23)  and  (28)  into  (32),  a  recurrence 
relation  for  F(AyB)  can  be  obtained.  For  t  €  Ay 

F(A,B)  =  PoF(A\t,B)  +  PjFiA(>B\j).  ■ 

/es\o 


The' algorithm  for  computing  F(A,B)  is  obtained  by  ^ 
recursive  implementation  of  (33).  The  initial  conditions 
associated  with  (33)  are  given  below.  '  '  .  •  ’ 


^  y  The  trajectory  of  target  i  i 
X  The  Jth  measurement  j 

<0  The  validation  gate  t 


Fig.  1.  Typical  layout  of  4  targets  and  4  measuremenis. 

1)  If  r  is  the  only  element  in  Ay  then, 

f(A,B)^'£Pj-  (3^ 

yes 

Note  that  Pq  =  Poin  the  above  equation  according  to 
the  notation  introduced  in  (11). 

2)  If  ;■  =  0  is  the  only  element  in  B,  then, 

F(A,B)  =  (Po)'^'  (3; 

where  |^|  denotes  the  number  of  elements  in  A. 

3)  If  1B|  =  2  and  \A\  >  \B\,  then, 

F(A,B)  =  (Po)'''!-'  {po+Yl  Pjl  •  (3f 


In  order  to  demonstrate  how  the  recursive 
algorithm  works,  consider  an  example. 

Example.  Suppose  that  there  are  4  closely  spaced 
targets  under  surveillance.  On  one  radar  scan,  four 
measurements  are  received  as  shown  in  Fig.  1.  Then, 

-  >1  =  {1,2,3, 4}  (3 

P  =  (0,1, 2, 3, 4}.  (3 

The  conventional  validation  matrix  occurring  in  the 
JPDAF  for  this  example  is  given  in  (39) 
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4. 

For  target  1,  there  is  only  one  measurement  which  falls 
inside  its  validation  gate.  Therefore,  only  /Jg  and 
are  none  zero.  Using  (23)  and  (28),  Pq  and  /3j  may  be 
computed  as 

/^^(A,B)  =  PoF(A\l,B)  (40) 

/3l(A,B)  =  P}F(A\l,B\l)  (41) 
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F({3,4},{0,2,3,4}) 


‘n{3,'l},{0,2,3}) 


^#•({4),  (0,2.3.41) 

^4  -f({4}.  (0.3,4)) 
F({4), {0.2.4}) 
F((4), (0.2.3)) 
^n(<).  {0.3.4}) 
F((l},(0.4}) 
F({4{.(0.3}) 


F({4),{0.2..1}) 
F({t},{0.4)) 
F()4),{0,2)) 
/;o,'F({4),  {0,2.3}) 


''/^({4).{0,2)) 

Fig.  2.  niustration  of  example  using  recursive  algorithm. 


where 

^\1  =  {2,3,4}  (42) 

B\1  =  {0,2,3,4}.  (43) 

Eqtiations  (40)  and  (41)  reduce  the  computations  of 
and  /?}  to  those  of  F(A\1,B)  and  F(A\1,B\1). 

Since,  in  this  example,  >1\1  =  {2,3,4}  and  5\1  = 
{0,2,3, 4},  the  computation  of  F(.4\1,B\1)  may  be 
performed  step  by  step  using  (33)  as  shown  through 
the  directed  tree  in  Fig.  2.  In  Rg.  2,  the  recursiw  • 
procedure  for  computing  F({2,3,4},  {0,2,3, 4})  is 
described  by  a  directed  tree.  In  this  tree,  the  Pjs  are 
the  weights  associated  with  the  edges.  The  plus  sign, 

+,  at  a  node  represents  the  summation  operation. 

The  F  function  at  a  node  labelled  with  a  plus  sign 
is  obtained  from  the  weighted  summation  of  the  F 
funaions  at  the  sons  (defined  to  be  nodes  which  have 
edges  incident  from  the  node  under  consideration)  of 
the  node,  as  seen  from  (33).  Let  a  leaf  denote  a  node 
which  does  not  have  sons.  Each  F  at  a  leaf  of  the  tree 
is  directly  computed  using  the  initial  condition  given 
in  (34).  In  the  computation  of  F({2,3,4},{0,^3,4}), 
in  this  example,-the  nodes  in  the  tree  shown  in  Rg.  2 
are  visited  via  DFS.  A  node  is  not  visited  if  the  only 
branch  incident  on  it  has  zero  weight  As  a  result,  the 
sons  of  this  node  are  also  not  visited.  In  Fig.  2,  the 
nodes  having  only  incident  dashed  branches  are  not 
visited.  In  other  words,  the  dashed  branches  in  the  tree 
are  pruned  becaiise'their  weights  are  zero. 

F(.4\1,F)  may.be  computed  in  a  similar  way  as 
discussed  above.  Generally  speaking,  the  computation 
in  the  recursive  algorithm  developed  here  is  performed 
from  top  to  bottom  (or  from  left  to  right  as  shown  in 
Fig.  2).  In  the  average  case  scenario,  a  lot  of  branches 
in  a  directed  tree  have  zero  weight,  such  as  the  one 


shown  in  Rg.  2.  The  computational  cost  can  be  saved 
if  a  directed  tree  is  visited  from  top  to  bottoni  since 
the  nodes  which  have  branches  with  zero  weights 
incident  onto  them  are  not  visited.  In  the  worst  case 
scenario,  however,  no  branches  have  zero  weight  ' 

The  overall  computation  may  be  reduced  if  the  nodes 
having  the  same  F  are  merged  so  that  they  are  visited 
only  once.  As  shown  in  Fig.  2,  almost  every  F  appears 
twice  at  the  leaves  of  the  tree  if  the  F  functions  at . 
the  nodes  which  have  branches  with  zero  weights 
incident  on  to  them  are  computed.  If  the  computation 
is  performed  from  bottom  to  top,  it  is  obvious  that 
the  computational  cost  can  be  reduced  in  the  worst 
case  scenario.  However,  in  the  average  case  scenario, 
the  majority  of  the  Fs  at  the  leaves  of  a  tree  appears 
only  once  after  the  dashed  branches  are  pruned,  as 
shown  in  Fig.  2.  Therefore,  a  lot  of  the  computation 
might  be  wasted  if  the  computation  is  performed  from 
bottom  to  top  since  many  nodes  which  are  connected 
with  dashed  branches  do  not  contribute  to  the  F  at 
the  root  of  a  tree.  It  is  shown  in  the  next  section  that 
the  computational  cost  is  lower  if  the  computation  is 
perform  from  top  to  bottom  in  the  average  case.  The 
computation  in  the  fast  JPDA  algorithm  [12]  may  be 
interpreted  as  being  performed  from  bottom  to  top  if 
the  procedure  is  organized  into  the  same  tree  structure 
as  discussed  above.  Detailed  comparison  between  the 
recursive  algorithm  developed  here  and  the  fast  JPDA 
algorithm  will  be  given  below. 


V.  COMPUTATIONAL  COMPLEXITY  ANALYSIS 

The  computational  complexity  measure  here  is 
defined  in  terms  of  the  numbers  of  multiplications, 
and  additions,  A(n,m),  used  for  calculating 
the  /3's  in  the  JPDAF,  where  n  is  the  number  of 
targets  and  m  is  the  number  of  measurements. 

Since  the  computational  cost  for  normalizing  /3y  is 
negligible  when  compared  with  that  of  F(A\t,B\j), 
multiplications  and  additions  for  normalizing  /3{  are 
not  included  in  and  A(n,fn).  In  the  worst 

case,  the  m  measurements  received  fall  inside  the 
intersection  of  the  validation  gates  of  the  n  targets. 
Therefore,  each  measurement  may  be  associated 
with  either  any  one  of  the  n  targets  or  clutter.  In  the 
average  case,  however,  it  is  assumed  that  there  are  on 
average  two  to  three  measurements  which  fall  inside 
the  validation  gate  of  each  target  In  the  following, 
the  worst  case  analysis  is  given  first  Then,  the  average 
case  analysis  is  discussed. 


A.  Worst  Case  Analysis 

In  the  worst  case,  no  branch  is  pruned  if  the 
computation  process  is  represented  in  a  tree  structure 
as  shown  in  Fig.  2,  without  any  dashed  edges. 

As  discussed  in  Section  III,  to  compute  all  P‘j, 
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only  F(A\t,B\j)  for  f  =  1,. and./ =  •• 

and  F(>t\l,5)  need  to  be  computed. .Therefore,  • 
M(n,m)  and  A(n,m)  are,  respectively,,  the  number 
of  multiplications  and  additions  which  suffice  for 
computing  nm  ^'jS  (j  jk  0)  and  one  In  the  recursive 
algorithm  proposed  in  the  previous  section,  the 
computation  of  the  jSjs  is  based  on  (23)  and  (28). 

As  shown  there,  the  only  operation  required  is 
one  multiplication  for  computing  each  /3y  after  the 
corresponding  F  is  computed.  Therefore, 

M(n,m)  =  nmM'(n  —  l,m) 

+M'(n-l,m  +  l)  +  nm  +  l  (44) 

and 

A{n,m)  =  nmA'(n  -l,m)+  A'(n  —  1,/n  + 1)  (45) 

where  M'(p,q)  and  A'(p,q)  denote,  respectively,  the 
numbers  of  multiplications  and  additions  which  suffice 
for  computing  the  function  F(A,B)  with  p  ■=\A\  and 
q  =  |5|. 

In  (33),  there  are  q  multiplications  and  q—1 
additions  used  in  the  computation  of  F{A,B)  after 
(q  —  V)  F(A\t,B\j)s  (J  ^  0)  and  one  F(A\t,B)  are 
computed.  Therefore, 

+  M'(p-l,q)’hq  (46) 

+  A'(p^l,q)  +  q^l  (47) 

with  boundary  conditions, 

M'(l,q)=0  .  (48) 

M'(p,2)  for  p>2  (49) 

M'(pA)^P-1  ..  -  ^  (50) 

A\hq)^q^l  ‘  J 

A'(p,2)-=p,  for  p>2  (52) 

^'(i7,l)=0.  ■  •  ^  v  -  (53) 

The  boundary  conditions  for  (46)  and  (47)  are 
obtained  from  (34),  (35),  and  (36). 

In  Thble  I,  some  sample  values  of  M(n,m)  ' 
andy4(n,m)  are  provided.  Comparing  the  listed 
values  in  Thble  I  with  the  corresponding  entries  in 
‘Tkble  II,  it  is  not  surprising  to  find  that  the  recursive 
algorithm  requires  more  multiplications  and  additions 
to  compute  the  /3yS  than  the  fast  JPDA  algorithm 
[12]  in  the  worst  case.  This  is  because  some  Fs  are 
computed  more  than  once  in  the  recursive  algorithm 
as  discussed  in  the  previous  section.  However,  the 
recursive  algorithm  is  expected  to  perform  much  better 
than  the  fast  JPDA  algorithm  in  the  average  case.  . 
rMore  discussion  on  the  average  case  analysis  will  be 
given  in  the  sequel. 


TABLE  I 

Sample  Values  of  M(n,m)/A(n,m)  in  the  Recursive  Algorithm 


11 

m 

iV/(n,  Tn)/A{n,  m) 

n 

m 

M{n,m)/A{7i,m) 

3 

6 

134/582 

5 

10 

43644/28390(1 

4 

2195/1  l!)r)2 

6 

12 

1034045/7807860 

TABLE  II 

Sample  Values  of  M(n,m)/A(n,m)  in  Fast  JPDA  Algorithm 


11 

_ 

m 

M{n,m)/ A{n^m) 

n 

in 

M(n,m)/A{7i.m) 

3 

6 

150/396 

5 

10 

10570/22100 

4 

8 

1404/3232 

6 

12 

70044/1365S4 

TABLE  III 

Sample  Values  of  M(n,m)/A(n,m)  in  DFS  Algorithm 


n 

m 

M(n,m)/A{n,m) 

n 

m 

M(n,m)/ 4(n,m) 

3 

6 

318/787 

5 

10 

96890/339041 

4 

8 

5072/14849 

6 

12 

2218992/9081085 

In  the  worst  case,  the  comparison  between  the 
computational  complexity  of  the  recursive  algorithm 
and  the  DFS  algorithm  is  also  based  on  M(n,m) 
and  A(nyni).  In  Table  III,  some  sample  values  of 
M{n,ni)  and  A(n,m)  required  for  implementing  the 
DFS  algorithm  are  listed.  Comparing  Thble  III  with 
Tkble  I,  it  can  be  inferred  that  the  M(n,m)s  in  the 
DFS  algorithm  are  more  than  twice  as  large  as  those 
in  the  recursive  algorithm  and  that  the  A(n,  m)s  in 
the  DFS  algorithm  are  always  larger  than  those  in 
the  recursive  algorithm.  Considering  the  fact  that  the 
recursive  algorithm  is  suitable  for  implementation  in  a 
multiprocessor  system,  the  computational  cost  of  the 
recursive  algorithm  could  be  much  less  than  that  of  the 
DFS  algorithm  in  the  worst  case. 

The  computational  complexity  given  above  is 
analyzed  in  terms  of  operational  counts.  It  can  also  be 
given  in  “big  O*"  notation.  If  the  algorithm  developed 
here  is  implemented  in  a  system  with  a  floating  point 
unit,  multiplication  and  addition  operations  would 
take  about  the  same  amount  of  time.  Let  G(/i,m)  = 
M(n,/n)  4-^(n,m).  According  to  the  definition  given 
in  [16,  p.  2]AVA,  the  computational  complexity  of  our 
algorithm  in  the  worst  case  scenario  is  0(G(n,m)). 

B.  Average  Case  Analysis 

In  order  to  simplify  the  discussion,  it  is  assumed 
that  there  are  at  most  three  measurements  inside  the 
validation  gate  of  each  target.  This  number  was  also 
selected  in  the  examples  given  in  [2,  8].  As  a  result, 
there  are  at  most  3n  /3js  (j  ^  0)  and  one  /3J  which 
need  be  computed.  After  the  interference  part  of  each 

is  obtained,  3n  4- 1  multiplications  are  sufficient 
for  computing  the  /3J.s  as  evident  from  (23)  and  (28). 
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table  rv  V 

Upper  Bounds  of  Af(n,m)  and  A(n,m)  in  the  Recursive  Algorithm 


a 

m 

Af(n,m)/i4(n,m) 

n 

m 

Af(n,m)/i4(n,m) 

.1 

- 

.50/150 

5 

- 

1360/4080 

4 

- 

273/819 

6 

~  i 

0A79119437 

Therefore,  (44)  and  (45)  lead  to  the  inequaUties, 
M(n,m)<3nM'{n-l,m) 

+  M'(n  -  l,m  + 1)  +  3n  + 1 


Upper  Bounds  of  Af  (n,m)M(rt,mym  DFS  Au 


n 

m 

Af(n,  m)/A(n,  m) 

q’ 

m 

3 

- 

90/208 

5 

1 

- 

417/1024 

6 

- 

5_  f  lL788/4^jf^:  : 


for  ready  reference  ' 

Af(rt,m)  <  2*4"  —  2- 3rt  —  3"  .(64) 

y4(rt,m)<4'’ +3n*4''-^  ’•-■*'■(65) 


A(n,m)  <  3nA'(n  -  l,m)  +  A'{n  -  \,m  +  1).  (55) 

In  the  average  case,  there  are  at  most  4 
multiplications  and  3  additions  required  in  the 
computation  of  F{A,B)  after  3  or  less  F(A\t,B\j')  s 
(J  ^  0)  and  one  F(A\t,B)  are  computed  in  (33). 

Again,  from  (46)  and  (47),  we  infer  that 

M'{p,q)  <  3M'(p  -  l,q  -  1)  •+•  M'(p  -  l,q)  +  4 

(56) 

A'(p,q)  <  2>A'{p  -  1,<?  - 1)  +^'CP  -  1.?)  +  3- 

(57) 

Since,  in  a  cluttered  enviromnent,  m  is  likely  to  be 
larger  than  n  in  the  average  case,  therefore  by  using 
(48)  and  (51),  the  initial  conditions  for  (56)  and  (57) 
may  be  simplified  to 

M'(1,?)  =  0  ,  (58) 


Comparing  the  upper  bounds  of  M{n,m)  in  (64) 
and  (62),  it  is  apparent  that  the  recursive  algorithm 
requires  less  multiplications  than  the  DFS  algorithm 
for  moderate  value  of  n.  A  similar  conclusion  may  be 
reached  by  comparing  the  upper  bounds  of  A{n,m)  in 
(63)  and  (65).  Furthermore,  we  infer  that  the  recursive 
algorithm  always  requires  less  additions  than  the 
DFS  algorithm.  The  computational  complexity  of  the 
recursive  algorithm  and  the  DFS  algorithm  in  the 
average  case  may  be  assessed  by  comparing  sample 
values  of  the  upper  bounds  of  M{n,m)  and  A(n,m) 
given  in  Table  IV  for  the  recursive  algorithm  and  in 
Table  V  for  the  DFS  algorithm. 

The  computational  complexity  of  our  algorithm  in 
the  average  case  can  also  be  given  in  terms  of  the  “big 
O”  notation.  From  (62)  and  (63),  one  can  conclude 
that  the  computational  complexity  of  our  algorithm  is 
0(n4''). 


A'(l,q)<3.  (59)  VI.  CONCLUDING  REMARKS 


With  the  above  initial  conditions,  it  is  not  difficult  to 
show  that  both  M'(p,q)  and  A'(p,q)  are  independent 
of  q  and  that  those  are  bounded  from  above  by  ■ 

M'(i7,9)<^(4^-*-l)  (60) 

A\p,q)<^‘’-1-  (61) 

Substituting  (60)  into  (54)  and  (61)  into  (55),  the  upper 
bounds  for  M{n,  m)  and  yl(n,  m)  in  the  average  case 
can  be  obtained  after  routine  manipulations.  _  ; 

M(/i,m)<5^(4'-i-l)  (62) 

A{n,m)  <  (3n  + 1)(4'-*  - 1).  (63) 

Since  the  design  of  the  fast  JPDA  algorithm 
[12]  is  based  on  the  worst  case  scenario,  there  is  no 
computational  cost  reduction  in  the  average  easel  As  a 
result,  the  recursive  algorithm  developed  here  is  much 
more  efficient  than  the  fast  JPDA  algorithm,  as  shown 
in  Thble  II  and  IV.  If  we  consider  the  fact  that  the 
entries  in  Tkble  IV  are  loose  upper  bounds,  the  actual 
computational  cost  of  the  recursive  algorithm  would  be 
much  less. 

The  upper  bounds  of  M(n,m)  and  in  the 

DFS  algorithm,  given  in  [13,  14],  are  reproduced  below 


In  this  paper,  a  recursive  algorithm,  which  is 
suitable  for  implementation  in  a  multiprocessor  system, 
is  developed.  In  this  algorithm,  the  computation  of 
the  a  posteriori  probabilities,  /3's,  is  not  based  on  the 
generation  of  the  data  association  hypotheses  like  in 
the  DFS  algorithm  [13].  Each  P]  in  this  algorithm  is 
decomposed  into  two  parts.  The  computation  of  one 
part  is  trivial  and  the  recursive  algorithm,  developed 
here,  is  used  to  compute  the  other  part  (due  to 
interference  from  other  targets). 

In  the  new  algorithm,  the  computation  of  the 
interference  part  of  is  implemented  recursively  in 
a  top-to-bottom  mode.  In  the  worst  case,  this  recursive 
algorithm  requires  more  multiplications  and  additions 
than  the  fast  JPDA  algorithm  [12].  However,  in  the 
average  case,  the  recursive  algorithm  is  expected  to 
outperform  the  fast  JPDA  algorithm.  In  comparison 
with  the  DFS  algorithm  developed  in  [13,  14]  the 
recursive  algorithm  requires  less  multiplications 
and  additions  than  the  DFS  algorithm.  The  most 
important  feature  of  the  recursive  algorithm  is  that 
it  can  be  implemented  on  a  multiprocessor  system. 
Some  suggestions  for  the  implementation  of  the  new 
algorithm  on  a  multiprocessor  system  are  given  in 
Appendix  A. 
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The  task  of  computing  F(A\t,B\j)  is  related  to 
that  of  evaluating  the  permanents,  defin^  j®  [17]>  of 
anmx  (n  —  1)  matrix  and  its  submatnces  having  Pf 
for  entries,  where  r  =  r  ^  f, 7  =  0,1,..., m,  and 

/  j.  Justifications  are  provided  in  Ajppendix  B,  where 
illustrative  examples  serve  to  verify  the  results  obtained 
by  exploiting  this  relationship. 

APPENDIX  A:  SUGGESTIONS' FOR  '  ^  ’ 
IMPLEMENTATION 

In  order  to  explore  the  inherited  parallelism  in 
our  efficient  algorithm,  the  computation  of  the  a 
posteriori  probabilities  in  our  algorithm  may 

be  summarized  as  below. 

1)  For  j  =  1,2,.. .,m,  t  =  l,2,...,n,  and  7^  0, 

^'j(A,B)  =  P‘jF(A\t,B\j)  (66) 

and  for  r  =  1,2,..., n, 

l3‘o  =  PoF(A\t,B).  (67) 

2)  For  any  A  and  B, 

F(A,B)  =  PoF(A\t,B)+ 

;6B\0 

(68) 

3)  The  constant 

c  =  ^P]{A,B)  -  (69) 

;=o 

may  be  computed  for  any  value  of  f  as  done  m  (29)  for 
the  r  =  l  ease.  ... 

4)  For  j  =  0,l,...,/n,  t  —  1,1,. ..,n,  and  ^ 

-  ^  '  '■■  ■'  ■  (70) 

To  implement  oiir  efficient  algorithm,  memory 
space  has  to  be  allocated  for  P‘j,  Pj,  A,  and  B.  Let 
BETA(0  :  m,  1 :  rt)  and  P(0  :  m,  1 :  n)  be  two  2-D 
arrays  each  of  size  (m  +  1)  x  n  for  storing  Pj  and 
Pj.  Since  A  and  B  are  ordered  sets  of  integers,  some 
special  data  structures  may  be  required  to  stored  these 
two  sets.  As  shown  in  Fig.  2,  two  types  of  operations 
are  required  on  A  and  5  to  compute  of  F(A,B). 

One  operator  deletes  an  element  from  either  Aor  B 
when  tihe  computation  advances  from  a' node  at  the 
higher  level  of  the  tree  to  a  node  at  the  lower  level 
The  other  operator  inserts  an  element  in  either  A 
or  B  when  the  computation  backtracks  from  a  node 
at  the  lower  level  of  the  tree  to  a  node  at  the  higher 
level.  Thus  each  operation  updates  either  A  or  B.  In 
order  to  implement  the  two  operations  efficiently,  data 
structures  can  be  found  in  [16]  for  storing  A  and  B. 

Let  the  two  functions,  Delete(.lf,;)  and  Insert(A',;), 


denote,  respectively,  the  operations  of  deletion  of 
element  j  from  set  X  and  insertion  of  element  j  into 
set  X  when  j  7^  0.  In  the  case  when  j  =  0,  Delete  and 
Insert  do  nothing.  Tb  compute  F(A,B),  a  function  can 
be  defined  using  (A.3).  For  the  sake  of  simplicity,  let 
F(A,B)  be  that  function  which  returns  the  value  on  the 
right-hand  side  of  (68). 

With  the  notations  introduced  above,  the 
computation  of  a  posteriori  probabilities  pj(A,B) 
in  our  efficient  algorithm  for  implementation  on  a 
multiprocessor  system  is  summarized  below. 

1)  For  ;  =0,1,2,... ,m,  t  =  1,2,...,«,  and  P(j,t) 

?^0, 

BETA(;,r)  =  P(;,r)F(Delete(A,t),Delete(5,;)). 

2)  The  normalization  constant  for  any  target  t  is 

m 

c  =  ^BETA(;,0. 

;=0 

3)  For  j  =  0,1,2,. t  =  1,2,...,«,  and  BETA(j.t) 

7^0, 

BETA(;,r)  iBETA(;,r). 

Each  BETA(/,r)  in  steps  1  and  3  above  can  be 
computed  in  parallel  in  a  multiprocessor  system 
such  as  CM-200  from  Connection  Machine  Inc.. 

A  programmer  may  view  the  CM-200  as  a  set  of 
virtual  processors^  one  for  each  data  element  The 
corresponding  code  and  data  are  passed  to  each 
processor.  On  a  CM-200,  steps  1  and  3  may  be 
implemented  using  a  conditional  FORALL  statement 
in  CM  FORTRAN  [18].  The  effect  of  a  conditional 
FORALL  statement  is  similar  to  the  embedding  of 
an  IF  statement  in  a  DO  loop  using  FORTRAN 
77.  With  a  conditional  FORALL  statement,  all  data 
elements  are  operated  on  simultaneously.  In  general, 
our  efficient  algorithm  is  ready  to  be  implemented 
in  any  parallel  computer  with  single-instruction  and 
multiple-data  (SIMD)  or  multiple-instruction  and 
multiple-data  (MIMD)  architecture. 

APPENDIX  B:  RELATIONSHIP  BETWEEN  F{A,B) 

AND  PERMANENT 

The  definition  for  a  permanent  is  given  below  [17]. 

Definition  Let  A  =  {afj)  be  an  m  x  n  matrix  over 
any  commutative  ring,  m  <n.  The  permanent  of  A, 
written  Per(^)  is  defined  by 

Per(A)  =  ^  (71  j 

<7 

where  the  summation  extends  over  each  of  the 
nl/(n  —  my,  one-to-one  maps,  denoted  by  a,  of  the  set 
{1,2,. ..,m}  to  the  set  (1,2,. ..,a2}. 
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In  our  case,  let  P  =  (Pj)  denote  the  matrix  havmg 
(m  + 1)  rows  and  n  columns.  A  and  B,  detoed  befor 
contain,  respectively,  the  column  and  row  indices 
of  P,  Le.,  A  =  {l»2,...,n}  and  B  =  {0,1,. .m'w}* 
notational  convenience,  in  later  usage,  we  denote  the 
permanent  of  P,  Per(P),  by  Per(>l,B).  This  causes 
no  confusion  or  ambiguity  because  the  2- tuple  (j>f) 
formed  from  an  element,  t,  of  A  and  an  element, 
of  B  uniquely  defines  Pj.  Therefore,  applying  the 
definition  for  permanent  to  a  matrix  or  its  transpose, 
we  get 

(T 

\B\<\A\  c 

?eT(A,  S)  =  K(ii)K(ii)  •  •  • 

<T 

\B\  >  Ml-  ( 

It  is  not  difficult  to  show  that  F(A,B)  in  (33)  is 
expandable  as 


equation 736). 


EXAMPLE  4  Let  ^  =  {i’3,4}  aha;B  ^ 

The  matrix  (P;),  where  the  indices^  andlf  are^orawn 

from  the  sets  ^  and  B,  is  (note  that  = 

/  p  p  p  \'  •  t  '  •• 
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Pi  Pi  Pi  ■  ;  ' 

PI  Pi  Pi  ■  ■ 

Pj  P|  Pt  :  ' 

VPj  p|  pV 

In  this  case,  \A\  <  |5\0|.  Then  (75)  applies,  yielding 
FiA,B)  =  (Po)^  +  (PofYl  E 

/6>t;6S\0 


\ieA 

+  (Po)'-^'“^[  Per({ti,«2}.S\0)  j 
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F(A,B)  can  also  be  calculated  as  below,  using  (33) 

F(A,B)  =  (PofJ2Pj+Po  E  E 

jeB  ;i€B\0 

+  '’oE^JEn 

h€B\0  /2G3\ii  V  ^  ^ 

+  En  Etf'E'>J- 

li6B\0  ;26B\{0,A}  ,  /56B\{/t.«} 


+  (Fo)l^'-^  f  E  .  ,  It  "is  not  difficult  to  show  that  (78)  and  (79)  are 

equivalent  .  * 

+  Per(A,0\O)  ,1S\0|  >1A1.  ■  ^  " 

The  vaUdities  of  (74)  and  (75)  are  iUustrated 

consideration  of  some  special  cases  below.  -  authors  tbanV  Jongtae  Chun  for  his  careful 

Example  1  Let  ^  =  {t}.  so  that  \A\  =  1.  In  this  case,  jea^ing  of  the  manuscript 
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A  Unified  Approach  to  Data  Association  in 
Multitarget  Tracking*t 

B.  ZHOUt  and  N.  K.  BOSE* 


Key  Words — Target  tracking;  tracking  systems. 


Abstract — A  unified  framework  is  proposed  to  provide  a 
comprehensive  understanding  of  the  problem  of  data 
association  in  multitarget  tracking.  Under  this  framework,  a 
systematic  scheme  is  developed  for  generating  the  data 
association  hypotheses  in  the  target-oriented,  measurement- 
oriented,  and  track-oriented  approaches.  Since  there  are 
many  data  association  hypotheses  with  identical  likelihoods 
in  the  measurement-oriented  and  track-oriented  approaches, 
two  specialized  algorithms  are  developed  to  efficiently 
generate  the  data  association  hypotheses  in  the  two 
approaches  by  adapting  the  depth-first  search  (DFS) 
algorithm  developed  earlier  in  the  target-oriented  case. 

1.  Introduction 

Recently,  considerable  attention  has  been  paid  to  the 
problem  of  data  association  in  multitarget  tracking.  In  a 
cluttered  environment,  the  detected  signals  comprising  the 
measurements  may  not  all  originate  from  the  targets  of 
interest.  Some  of  those  could  be  either  from  clutter  or  false 
alarm.  Each  detected  signal  may  be  used  as  a  measurement 
to  update  the  state  of  a  target.  Therefore,  more  than  one 
measurement  may  be  available  for  the  purpose  of  target  state 
updating.  The  ambiguity  occurring  in  associating  previously 
known  targets  with  new  measurements  is  referred  to  as  the 
problem  of  data  association.  To  solve  this  problem,  three 
approaches  have  been  reported  in  the  literature.  Those  are 
the  target-oriented,  measurement-oriented,  and  track -oriented 
approaches.  In  the  target-oriented  approach,  each  measure¬ 
ment  is  assumed  to  have  originated  from  either  a  known 
target  or  clutter  (Bar-Shalom  and  Fortmann,  1988; 
Bar-Shalom  and  Tse,  1975;  Fortmann  et  ai  1983).  In  the 
measurement-oriented  approach,  each  measurement  is 
hypothesized  to  have  originated  from  either  a  known  target, 
a  new  target,  or  clutter  (Reid,  1979).  In  the  track -oriented 
approach,  however,  each  track  is  hypothesized  to  be  either 
undetected,  terminated,  associated  with  a  measurement,  or 
linked  to  the  start  of  a  maneuver  (Kurien  and  Washburn, 
1985;  Kurien  and  Liggins,  1988;  Bar-Shalom,  1990).  In 
multitarget  tracking,  Fortmann  et  ai  (1983)  proposed  the 
Joint  Probabilistic  Data  Association  Filter  (JPDAF)  as  an 
implementation  of  the  target-oriented  approach.  Since  the 
JPDAF,  which  may  be  viewed  to  be  synonymous  with  the 
target-oriented  approach,  lacks  track  initiation  capability  and 
also  suffers  from  the  track  biases  and  coalescence  problem 
(Bar-Shalom,  1990;  Fitzgerald,  1986),  further  study  of  the 
measurement-oriented  and  the  track-oriented  approaches 
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becomes  necessary  to  address  the  computation  intensive 
algorithms  associated  with  these  approaches  in  comparison  to 
the  JPDAF.  Let  n  and  m  denote,  respectively,  the  number  of 
targets  and  the  number  of  measurements.  In  the  worst  case, 
each  of  the  n  targets  could  be  associated  with  any  one  of  the 
m  measurements.  It  has  been  shown  that  the  computational 
cost  for  associating  the  n  targets  with  the  m  measurements 
increases  exponentially  as  a  function  of  the  parameters  n  and 
m  (Zhou  and  Bose,  1993).  However,  if  the  data  association 
problem  involving  n  targets  and  m  measurements  can  be 
decomposed  into  smaller  subproblems,  then  the  overall 
computational  cost  may,  approximately,  become  polyno- 
mially  dependent  on  n  and  m.  The  resulting  small  groups  of 
targets  and  measurements  are  referred  to  as  clusters  in  the 
multitarget  tracking  terminology.  Two  targets  are  in  the  same 
cluster  if  and  only  if  at  least  one  measurement  falls  inside  the 
intersection  of  the  validation  gates  of  the  two  targets.  A 
validation  gate  of  a  target  is  a  region  of  certain  size  that 
surrounds  the  predicted  position  of  the  target.  In  this  paper, 
the  data  association  problem  is  with  reference  to 
measurements  and  targets  belonging  to  the  same  cluster. 

After  targets  are  properly  clustered,  Zhou  and  Bose  (1993) 
proposed  a  depth-first  search  (DFS)  algorithm  to  efficiently 
generate  the  data  association  hypotheses  in  the  JPDAF. 
They  showed  that  a  significant  reduction  in  the  computa¬ 
tional  cost  of  implementation  of  the  JPDAF  can  be  realized 
by  the  sharing  of  common  factors  occurring  in  the  calculation 
of  the  conditional  probabilities  of  the  data  association 
hypotheses.  In  this  paper,  the  problems  of  generation  of  data 
association  hypotheses  in  the  measurement-oriented  and 
track-oriented  approaches  are  reformulated  to  permit 
development  of  a  mathematical  model,  similar  to  the  one 
developed  by  Zhou  and  Bose  (1993)  for  the  data  association 
problems  in  the  target-oriented  case.  A  unified  framework 
is  provided  for  adapting  the  DFS  algorithm  to  efficiently 
generate  the  data  association  hypotheses  in  the 
measurement-oriented  and  track -oriented  approaches,  where 
the  total  numbers  of  data  association  hypotheses  are 
expected  to  increase  drastically  over  the  target-oriented 
approach.  However,  reduction  in  the  overall  computational 
cost  may  be  feasible  from  observations  on  the  conditional 
likelihood  of  a  data  association  hypothesis.  In  the 
target-oriented  approach,  the  conditional  likelihood  of  each 
data  association  hypothesis  is  unique.  When  targets  are 
grouped  into  clusters,  this  uniqueness  property  does  not  hold 
for  the  measurement-oriented  and  track-oriented  ap¬ 
proaches.  Our  goal  here  is  to  develop  two  specialized  DFS 
algorithms  which  can  efficiently  identify  the  data  association 
hypotheses  with  identical  conditional  likelihood  in  the 
measurement-oriented  and  track -oriented  approaches. 


2.  Problem  formulation 

Suppose  that  there  are  two  targets  moving  towards  each 
other.  In  one  radar  scan,  three  measurements  are  received, 
as  shown  in  Fig.  1.  In  the  target-oriented  approach,  a 
measurement  may  originate  from  either  a  known  target  or 
clutter.  Zhou  and  Bose  (1993)  described  this  situation  using 
three  finite  sets  of  ordered  integers,  (;  -  L  2.  3).  The  three 
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sets  are 

2i={0,  1},  Z2  =  {0,1,2},  Z3  =  {0.2}. 

Zi  =  (0,  1}  implies  that  measurement  1  could  originate  from 
either  clutter  or  target  1.  Similar  explanations  may  be  given 
to  Z2  and  Z3.  In  the  measurement-oriented  approach, 
however,  any  one  of  the  three  measurements  could  have 
originated  from  a  new  target.  The  three-finite  sets  of  ordered 
integers  may  be  augmented  in  the  following  manner  to 
incorporate  the  possibility  that  a  measurement  may  originate 
from  a  new  target: 

Z,  =  (0,  1,  3},  Z2  -  (0, 1,  2,  4},  Z3  =  (0,  2,  5}. 

Note  that  Z,- {0,1,3}  implies  that  measurement  1  could 
originate  either  from  clutter,  or  old  target  1,  or  a  new  target 
3.  Similar  comments  apply  to  Z2  and  Z3.  In  addition  to  the 
possibility  that  a  measurement  could  originate  from  a  new 
target,  an  established  track  could  either  be  terminated  or 
linked  to  the  start  of  a  maneuver  in  the  track-oriented 
approach.  For  the  sake  of  simplicity,  suppose  that  Singer’s 
djmamic  model  (Singer,  1970)  is  used.  Since  the  man¬ 
euverability  of  a  target  is  inherited  in  Singer’s  dynamic 
model,  it  is  not  necessary  to  consider  it  in  data  association. 
On  the  other  hand,  a  new  integer  set  can  be  created  to  allow 
for  the  possibility  that  a  track  may  be  terminated  at  any 
moment.  For  the  example  shown  in  Fig.  1,  four  sets  of  finite 
ordered  integers  may  be  used: 

Zo  =  {l,2},  Zi={0,l,3},  Z2-{0, 1,2,4},  Z3  =  {0,2,5}, 

where  Zo-{l,2}  implies  that  the  tracks  of  targets  1  and  2 
may  be  terminated.  Note  that,  as  in  the  measurement- 
oriented  approach,  Z2  =  {0, 1,  2,  4}  implies  that  measurement 
2  could  originate  either  from  clutter,  old  target  1  or  2,  or  a 
new  target  4. 

By  comparing  the  data  association  problems  in  the 
measurement-oriented  and  track-oriented  approaches  with 
their  counterpart  in  the  JPDAF,  the  DFS  algorithm 
developed  by  Zhou  and  Bose  (1993)  may  be  easily  extended 
to  generate  the  data  association  hypotheses  in  the 
measurement  and  track-oriented  approaches.  Since  a 
measurement  or  a  track  may  be  associated  with  more  choices 
in  these  two  approaches  in  contrast  to  the  target-oriented 
case,  the  total  number  of  data  association  hypotheses  is 
expected  to  increase.  Therefore,  it  is  important  to  reduce  the 
computational  cost  in  such  situations. 

In  the  target-oriented  approach,  the  conditional  likelihood 
of  each  data  association  hypothesis  is  unique.  This 
uniqueness  property  does  not  hold,  in  general,  in  the 
measurenient-oriented  and  track-oriented  approaches.  Con¬ 
sider,  again,  the  example  described  in  Fig.  1.  Suppose  that 
measurements  1,  2  and  3  originate,  respectively,  from  target 
1,  clutter,  and  a  new  target.  Then,  the  conditional  likelihood 
of  this  data  association  hypothesis  is  (Reid,  1979): 

(1) 

where  c  is  a  constant,  F, ,  is  the  likelihood  that  measurement 
1  originates  from  target  1,  is  the  density  of  clutter,  and 


A^t-  is  the  density  of  previously  unknown  targets. 
Alternatively,  measurement  2  could  originate  from  a  new 
target  and  measurement  3  could  originate  from  clutter.  It  can 
be  verified  that  the  conditional  likelihood  of  this  data 
association  hypothesis  equals  the  one  given  in  (1).  Therefore, 
there  are  at  least  two  data  association  hypotheses  with  equal 
conditional  likelihood  in  this  example.  In  general,  A^  and 
A^t-  may  not  be  constant  in  a  surveillance  region.  For 
example,  clutter  distribution  may  not  be  uniform  and  the 
measurements  which  lie  in  the  center  of  the  surveillance 
region  (or  field-of-view)  are  less  likely  to  come  from  new 
targets  than  those  which  lie  near  the  boundary  of  the 
surveillance  region.  However,  measurements  and  targets  are 
normally  grouped  into  clusters  in  multitarget  tracking. 
Within  each  cluster,  A^  and  A  vt-  can  be  treated  as  constants. 
When  the  number  of  targets  and  the  number  of 
measurements  in  a  cluster  are  large,  duplications  in 
conditional  likelihood  can  occur  much  more  often.  The  same 
observation  can  be  made  in  the  track-oriented  approach. 
Therefore,  specialized  DFS  algorithms  which  can  effectively 
identify  such  duplications  in  the  measurement-oriented  and 
track-oriented  approaches  need  to  be  developed  to  reduce 
computational  cost  significantly. 


3.  Development  of  an  unified  framework 

In  the  target-oriented  approach,  the  received  measure¬ 
ments  are  divided  into  two  groups.  The  measurements  in  the 
first  group  are  assumed  to  be  associated  with  some  known 
targets  and  the  measurements  in  the  second  group  are 
assumed  to  be  from  clutter.  The  data  association  between  the 
known  targets  and  the  measurements  in  the  first  group  is 
accomplished  by  the  DFS  algorithm  (Zhou  and  Bose,  1993). 
The  association  between  clutter  and  the  measurements  in  the 
second  group  is  trivial. 

In  the  measurement-oriented  approach,  a  measurement 
could  be  associated  with  either  a  known  target  or  clutter  or  a 
new  target.  To  accommodate  these  three  possibilities,  the 
received  measurements  may  be  divided  into  three  groups.  As 
pointed  out  in  the  previous  section,  some  data  association 
hypotheses  in  the  measurement-oriented  approach  may  have 
identical  conditional  likelihood.  If  those  hypotheses  can  be 
identified,  the  corresponding  likelihood  can  then  be 
computed  only  once.  Suppose  that  there  are  m 
measurements.  In  a  particular  data  association  hypothesis  e, 
measurements  are  associated  with  previously  known 
targets,  measurements  are  associated  with  clutter,  and 
N^r  measurements  are  associated  with  new  targets. 
Therefore,  m  —  +  Nc  +  N/vt-.  It  can  be  shown  that  a  factor 

of  the  conditional  likelihood  of  e  which  corresponds  to  the 
association  between  measurements  and  clutter  is 

proportional  to  A^*^.  Similarly,  a  factor  of  the  conditional 
likelihood  of  e  which  corresponds  to  the  association  between 
measurements  and  new  targets  is  proportional  to 
When  a  new  data  association  hypothesis  is  generated  by 
merely  shuffling  the  measurements  associated  with  clutter  in 
£  with  those  associated  with  new  targets,  the  conditional 
likelihood  of  this  new  hypothesis  will  be  the  same  as  that  of  e 
as  long  as  Nc  and  are  unchanged.  This  observation 
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suggests  that  the  received  measurements  should  be 
partitioned  in  a  certain  way  for  the  purpose  of  identifying 
data  association  hypotheses  with  identical  conditional 
likelihood. 

The  difference  between  the  track-oriented  and  measurement- 
oriented  approaches  is  in  the  manner  of  handling  of  the 
previously  known  targets  which  are  not  associated  with  any 
measurements  in  a  data  association  hypothesis.  In  the 
measurement-oriented  approach,  such  targets  are  simply 
classified  as  undetected.  However,  in  the  track -oriented 
approach,  those  targets  could  be  either  undetected  or  their 
trajectories  could  have  been  terminated.  In  other  words,  the 
previously  known  targets  are  divided  into  three  distinct 
groups  in  the  track -oriented  approach.  The  partitioning  of 
the  previously  known  targets  can  be  done  in  the  same  way  as 
the  partitioning  of  the  measurements  in  the  measurement- 
oriented  approach. 

In  conclusion,  the  generation  of  the  data  association 
hypotheses  in  all  three  approaches  may  be  described  as 
below. 

1.  Select  the  measurements  and  the  previously  known  targets 
which  are  to  be  associated  with  each  other.  If  the 
target-oriented  approach  is  used,  then  a  data  association 
hypothesis  is  generated  and  this  process  is  repeated  until 
all  hypotheses  are  generated.  Otherwise,  go  to  step  2. 

2.  Partition  the  remaining  measurements  into  two  groups. 
The  measurements  in  one  group  are  hypothesized  to  be 
associated  with  clutter.  The  measurements  in  the  other 
one  are  then  to  be  associated  with  new  targets.  If  the 
measurement-oriented  approach  is  used,  then  go  back  to 
step  1  and  this  process  is  repeated  until  all  data 
association  hypotheses  are  generated.  Otherwise,  go  to 
step  3. 

3.  Partition  the  remaining  known  targets  into  two  groups. 
The  targets  in  one  group  are  hypothesized  to  be 
undetected.  The  targets  in  the  other  group  have  their 
trajectories  terminated.  Repeat  steps  1-3  until  all  data 
association  hypotheses  are  generated. 

The  above  description  provides  a  unified  framework  for  the 
generation  of  the  data  association  hypotheses  in  all  three 
approaches.  It  can  be  claimed  that  the  measurement-oriented 
approach  is  a  special  case  of  the  track-oriented  approach 
when  no  target  has  its  trajectory  terminated  and  the 
target-oriented  approach  is  a  special  case  of  the 
measurement-oriented  approach  when  there  are  no  new 
targets.  Under  this  unified  framework,  the  association 
between  the  measurements  and  the  previously  known  targets 
can  be  established  by  adapting  the  DFS  algorithm  originally 
developed  for  the  target-oriented  case  (Zhou  and  Bose, 
1993). 

4.  Two  DFS  algorithms 

The  data  association  problem  in  the  measurement-oriented 
approach,  as  discussed  in  the  last  section,  can  be  solved  in 
two  stages.  First,  the  association  between  the  measurements 
and  the  previously  known  targets  can  be  established  using 
the  DFS  algorithm  (Zhou  and  Bose,  1993).  Then,  the 
remaining  measurements  are  divided  into  two  groups.  The 
measurements  in  one  group  are  assumed  to  be  associated 
with  clutter,  and  those  in  the  other  group  are  assumed  to 
originate  from  new  targets.  Nijenhuis  and  Wilf  (1975) 
presented  a  partitioning  algorithm  which  can  efficiently  draw 
different  samples  of  size  k  out  of  N  objects.  After  Nr 
measurements  are  associated  with  Nj  targets,  Nc 
measurements  may  be  selected  from  the  (m  Nr)  remaining 
measurements  for  association  with  clutter.  Subsequently, 
Nf^r  ^(m-  Nr-  Nc)  measurements  are  assumed  to  be  from 
new  targets.  By  cascading  the  DFS  algorithms  and  the 
partitioning  algorithm,  a  new  DFS  algorithm  may  be 
developed  to  generate  all  data  association  hypotheses  in  the 
measurement-oriented  approach.  After  each  association 
between  Nr  measurements  and  Nr  targets  is  established, 
for  a  given  Nc  {  —  0, . . , ,  m  —  A/7 ),  the  likelihood 
Ajjf is  computed  only  once.  However,  for  a 
given  Nc,  (m  -  NrVJiNcH^  ^  -  ^<  )^>)  different  parti¬ 


tions  can  be  generated,  where  K\  denotes  the  factorial  of  K. 
This  means  that  the  (m  -  A/7  )!/(A/c-!(m  -  A/^  -  A/^^  )!)  data 
association  hypotheses  with  identical  likelihood  are 
effectively  grouped  together  in  this  new  DFS  algorithm. 

In  the  track-oriented  approach,  both  the  remaining 
measurements  and  the  remaining  targets  are  divided  into  two 
groups  after  Nr  measurements  are  associated  with  A/7-  targets. 
Suppose  that  A/^,  targets  have  their  tracks  terminated  (or 
discontinued)  and  A/^vd  targets  are  not  detected.  If  there  are 
n  previously  known  targets,  then  n=A/7-  +  A/o  +  A/^£)- 
Similar  to  the  measurement-oriented  approach,  a  new  DFS 
algorithm  may  be  developed  by  cascading  the  DFS  algorithm 
in  the  target-oriented  case  (Zhou  and  Bose,  1993)  and  two 
partitioning  algorithms  (Nijenhuis  and  Wilf,  1975).  In  the 
track-oriented  approach,  there  could  be 

_ {m-Nr)\{n-Nr)\ _ 

Ncl  (m-Nr-  A/c)!  N^l  (m-Nr-  No)\ 

data  approach  hypotheses  with  identical  likelihood. 

As  shown  above,  a  lot  of  computation  may  be  saved 
by  using  the  two  new  DFS  algorithms  to  generate  data 
association  hypotheses  in  the  measurement-oriented  and 
track -oriented  approaches.  The  computational  complexity  of 
the  two  algorithms  may  be  analyzed  as  was  done  for  the  DFS 
algorithm  in  the  target-oriented  case  (Zhou  and  Bose,  1993). 

5.  Conclusion  and  remarks 

The  unified  framework  developed  in  this  paper  provides  a 
systematic  scheme  for  the  generation  of  data  association 
hypotheses  in  the  target-oriented,  measurement-oriented  and 
track-oriented  approaches.  Furthermore,  the  framework 
makes  the  development  of  efficient  algorithms  possible  for 
the  measurement-oriented  and  track-oriented  approaches  by 
adapting  the  DFS  algorithm  developed  in  the  target-oriented 
scenario.  The  success  of  the  algorithms,  in  the  measurement- 
oriented  and  track-oriented  approaches,  is  derived  from 
exploitation  of  the  fact  that,  after  the  targets  and 
measurements  are  grouped  into  clusters,  many  data 
association  hypotheses  share  an  identical  conditional 
likelihood.  The  computational  cost  can  be  saved  significantly 
by  computing  the  identical  likelihood  only  once  in  each  such 
situation. 

The  ambiguity  of  data  association  in  the  measurement- 
oriented  and  track-oriented  approaches  may  be  further 
resolved  using  the  data  in  the  N  subsequent  scans.  Singer  et 
ai  (1974)  showed  that  near-optimal  performance  was 
achieved  with  N=l  for  the  single  target  case.  In  the  case  of 
multitarget  tracking,  the  data  from  two  successive  scans 
{N-2)  were  used  in  data  association  (Bar-Shalom,  1990). 
The  process  of  using  data  from  more  than  one  scan  in  data 
association  is  called  multiscan  correlation.  The  application  of 
the  algorithms  developed  in  this  paper  to  multiscan 
correlation  is  straightforward.  To  generate  data  association 
hypotheses  in  multiscan  correlation  in  the  measurement- 
oriented  and  track-oriented  approaches,  the  algorithms 
proposed  in  Section  4  may  be  applied  to  generate  data 
association  hypotheses  in  each  scan. 
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niques  for  parameter  estimation  of  multiple  signals  in  addi¬ 
tive  white  noise  or  the  directions  of  arrivals  of  wavefronts 
at  an  array  of  sensors  use  either  the  information  generated 
from  the  noise  subspace  or  the  signal  subspace.  Here,  we 
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other  subspace  techniques  are  made  and  an  illustrative  ex¬ 
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Since  1979,  several  subspace  methods  for  parame¬ 
ter  estimation  have  been  introduced.  MUSIC,  ES¬ 
PRIT,  and  their  variants  such  as  TLS-ESPRIT  and 
GEESE,  which  can  be  applied  to  estimate  the  direc¬ 
tions  of  arrivals  of  impinging  wavefronts  by  operating 
on  the  signals  collected  at  an  array  of  sensors,  are 
described  in  [1] .  More  recently,  the  weighted  sub¬ 
space  method  introduced  in  [2]  has  been  shown  to  be 
asymptotically  efficient  on  the  estimation  error  vari¬ 
ance,  as  the  stochastic  maximum  likelihood  technique 
is  known  to  be,  under  certain  regularity  conditions 
involving  the  assumption  of  Gaussian-distributed 
emitter  signals. 

Though  subspace  methods  provide  only  one  para¬ 
digm  for  signal  processing,  their  popularity  and  suc¬ 
cess  in  implementation  may  be  attributed  to  the  sim- 
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pie  underlying  model  for  the  signal  and  noise.  The 
assumptions  on  the  model  could  range  from  the  re¬ 
quirement  of  having  a  smaller  number  of  signals  than 
sensors  to  a  noise  field  characterization  which  may  be 
unrealistic  for  real  data  gathered  from  sensors.  Never¬ 
theless,  the  subspace  methods  provide  a  simple  uni¬ 
fied  approach  to  tackle  inherently  difficult  nonlinear 
parameter  estimation  problems,  and  their  scopes  con¬ 
tinue  to  increase.  For  a  recent  review  of  approaches  to 
tackle  the  nonlinear  optimization  problem  which 
yield  array  localization  signal  processing  techniques 
for  estimating  the  number  as  well  as  parameters  of 
narrowband  and  wideband  sources,  see  [3] .  For  a  re¬ 
cent  review,  in  a  unifying  subspace -fitting  framework 
of  methodologies  for  analyzing  a  variety  of  subspace 
algorithms,  see  [  4  ] . 

In  this  paper,  we  start  with  a  model  of  second-order 
data  (correlation  matrix)  generated  from  the  samples 
of  P  received  signals  { Xj  (n) } ,  i  =  1,  2, . .  . ,  P.  Each  of 
these  signals  contains  either  a  known  or  estimated 
number  M  of  moderately  correlated  signals  ( n ) ,  = 
1,  2,  ...  ,  My  having  distinct  directions  of  arrivals  6^, 
02,  .  .  •  ,  with  respect  to  a  linear  array  of  point  sen¬ 
sors,  embedded  in  additive  independent  and  identi¬ 
cally  distributed  Gaussian  noise,  as  shown  in  ( 1 )  be¬ 
low.  The  assumption  that  the  signal  subspace  dimen¬ 
sion  is  either  known  or  estimated  is  made  in  all 
subspace -based  methods  including  MUSIC,  ESPRIT, 
and  GEESE.  The  dimension  of  the  signal  subspace  is 
estimated  by  comparing  the  magnitudes  of  the  eigen¬ 
values  of  the  correlation  inatrix  and  using  a  threshold. 
For  i  -  1,  2,  .  .  .  ,  P,  7  =  V^,  and  =  tt  cos  0^, 

M 

Xi(n)  =  Z  3'fe(^)exp[-yco*(i  -  1)]  +  Vi(n), 
k^l 

n=  1,2,  (1) 
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It  is  assumed  that 


the  matrix  A  can  be  rewritten  as 


(2) 

and  noises  at  each  of  the  P  sensors  are,  ideally,  zero- 
mean,  wide- sense  stationary  uncorrelated  random 
processes.  Suppose  that  the  common  variance  of  the 
identically  distributed  uncorrelated  noise  processes  is 
Define 

x(n)  =  [x^(n)  X2(n)  •  •  •  x^in)]^  (3a) 

and 

y(n)  =  lyi(n)  y2(n)  •  •  •  yM(n)]\  (3b) 

where  the  superscript  “t”  denotes  “transpose.”  The 
standard  model  for  the  (P  X  P)  correlation  matrix 
[1], 

S  =  B[x(n)x^(n)],  (4) 


A  =  [q{e  ''“^)a(c  ■'“2)  .  .  .  a(e  ''^)]. 

Since  the  co/s  for  i  =  1,  2, .  .  . ,  Mare  distinct,  the  {P  X 
M)  matrix  A  is  of  full  rank  M,  and  because  P  is  a 
nonsingular  {M  X  M)  Hermitian  matrix,  the  matrix 
ARA^  must  be  a  (P  X  P)  nonnegative  definite  Her¬ 
mitian  matrix  of  rank  M.  Therefore,  ARA^  has  M 
nonzero  real  eigenvalues  /z, ,  f  =  1, 2, .  .  .  ,  M.  Since  S  is 
Hermitian  and  is  modeled  as  in  ( 5 ) ,  it  must  be  diag- 
onalizable  by  a  unitary  matrix  B  [5].  Therefore, 

B^SB  =  B^{ARA^-\-  a^I)B 

=  Diag[(Ati  +  cr^)  •  •  •  (mm+ *  •  <^]  (8) 

is  a  (P  X  P)  diagonal  matrix.  For  notational  brevity, 
we  define  below  the  (M  X  M)  and  (P  —  M)  X  (P  —  M) 
diagonal  matrices,  Ai  and  Ag,  respectively: 


where  the  superscript  “H”  denotes  “complex  conju¬ 
gate  transpose,”  is 

S  =  ARA^  (5) 

which  is  a  (P  X  P)  matrix,  R  =  E[y{n)  y^{n)]  is  an 
(M  X  M)  nonsingular  Hermitian  matrix,  and 


A  = 


1 


e 


1 


(6) 


The  matrix  I  denotes  throughout  an  identity  matrix 
of  appropriate  order;  to  avoid  notational  clutter,  we  do 
not  show  the  order  explicitly  in  /.  The  above  model 
provided  in  (1)  with  respect  to  the  direction -of-ar- 
rival  ( DOA )  problem  applies  also  in  several  other  situ¬ 
ations,  including  the  problem  of  parameter  estima¬ 
tion  of  sinusoids. 


Ai  =  Diag[  XjXz  •  •  •  Xj(^] ,  X;  =  m,-  + 

for  i  =  (9a) 

^2  ~  Diag[  Xa^+iXjv^+2  '  ■  ■  ^p]  >  ^M+i  “ 

for  i=l,...,P-M.  (9b) 

Denote  the  (P  X  M)  matrix  composed  of  the  eigen¬ 
vectors  associated  with  the  eigenvalues  Xi ,  Xg , . .  . , 
by  Ps  (signal  subspace  eigenvector  matrix)  and  the 
(P  X  (P  —  M))  matrix  composed  of  the  (P  —  M) 
eigenvectors  associated  with  the  eigenvalues  X^ ,  /  =  M 
+  1,  M  +  2, .  . .  ,  P  by  Pn  (noise  subspace  eigenvector 
matrix ) ,  so  that 


P  =  [PsPnJ.  (10) 

Some  well-known  results  are  rederived  to  facilitate 
the  documentation  of  the  novel  results.  From  (10) 
and  (8),  since  P”P  =  7, 


2.  MAIN  RESULTS 

2.7.  Analysis  of  the  Standard  Model 

The  standard  model  for  the  (P  X  P)  correlation 
matrix  S  was  given  in  (5).  From  (4)  and  (5), 

S  =  E[xin)x^{n)]  =  ARA^  -h  a-"/, 

where  A  is  shown  in  ( 6 ) .  On  defining 

q{e~-^^^)  =  [1  (7) 


[BsB^]^(ARA^)[BsB^]  +  =  Diag[A,  A^].  (11) 

The  above  equation  leads  to  the  following  two  equali¬ 
ties: 

APA^Ps  =  PgAi  -  B^ia^I) 

=  PsDiag[Mi*  •  -MAf]  (12a) 

APA^Pj^j  =  PJ4A2  Pjvj(o’^7)  =  0.  (12b) 

Obviously,  in  ( 12a)  and  ( 12b) ,  the  /’s  denote  identity 
matrices  of  orders  M  and  (P  -  M),  respectively. 


Since  AR  is  of  rank  M,  it  follows  from  ( 12b)  that  the 
(MX  (P  -  M))  matrix  must  satisfy 

(13) 

Since  Ps^s  (12a)  implies 

(Ps  A)P(A^Ps)  =  Diag[iUi,  ^^2^  •  •  •  »  Mm]  » 

which  in  turn  implies  that  the  (M  X  M)  matrix  A  ^Ps 
is  of  nonsingular.  Thus,  it  follows  from  (12a)  that 

A  =PsDiag[Mi-  •  =  BsC.  (14) 

where  C  =  Diag[)Lti  •  •  •  MmK^^^s)  ^  obvi¬ 
ously,  nonsingular.  It  is  seen  later  that  the  columns  of 
C  are  related  to  the  eigenvectors  of  the  matrix  whose 
eigenvalues  estimate  the  parameters  in  the  standard 
model  of  the  signal  correlation  matrix. 

Alternatives  to  ( 13 )  and  (14)  may  also  be  obtained 
from  (5) .  Since  in  practice  the  matrix  P  violates  the 
assumptions  of  the  data  model,  the  use  of  a  modified 
P  ( like  the  matrix  obtained  from  P  by  complex  conju¬ 
gation  of  its  elements)  is  likely  to  amplify  the  errors 
in  parameter  estimation.  Therefore,  the  alternate 
derivation  given  next  is  necessary,  where  the  counter¬ 
parts  of  (13)  and  (14)  are  obtained  in  (20)  and  (21) 
below  directly  from  the  model.  This  type  of  derivation 
minimizes  the  ambiguity  that  is  likely  for  the  location 
of  the  roots  of  a  polynomial  to  be  on  the  unit  circle 
when  certain  necessary  conditions  on  the  coefficients 
of  the  polynomial  are  not  satisfied  as  explained  later 
in  Subsection  2.2  below.  Let  the  exchange  matrix  be 
denoted  by 

’  0  0 
0  0 

:  : 

6  1 

1  0 

Premultiplying  and  postmultip lying  both  sides  of  (5) 
by  J,  we  have 

JSJ  =  E[{Jx{n)}  {Jx{n)}^] 

=  (JA)R{JA)^  +  a^I.  (15) 

JA  is  rearranged  as 

JA  =  •  •  Ja(e‘^“")] 

=  [a(e-'"')-  •  •a(e-'“’")] 


X  Diag[e”-'‘^‘^’“‘  •  •  •  e 
=  A+A, 

where 

a(e-'“')  =  [1  * 

A+  =  [a(e''“')-  •  •a(e^"")], 

and 

A  =  Diag[g"-'‘^“”“'  •  •  • 

Using  the  above  equation,  (15)  is  rewritten  as 

JSJ  =  A+(AfiA”)A?  +  =  A^R^A'^  +  <r"/,  (16) 

where  R+  =  Ai?A“  is,  obviously,  nonsingular. 

The  matrix  JSJ  is  Hermitian  and  since 

B^SB  =  (JB)^(JSJ}{JB) 

=  Diag[XjX2'  •  ■  ’^p]  (1"^) 

with  (JB)^{JB)  =  (JB)(JB)^  =  I,  therefore  JB  is 
the  eigenvector  matrix  associated  with  eigendecom- 
position  of  JSJ.  Furthermore,  since  JB  =  [JBsJBj^] 
=  [ Bsjnj]  where  Bgj  =  JBq  and  Bj^j  =  JBn*  the  (P  X 
M)  eigenvector  matrix  Bgj  is  associated  with  X/s,  i  = 
1,  2,  .  .  .  ,  M  and  the  (P  X  (P  -  M) )  eigenvector  ma¬ 
trix  Bnj  is  associated  with  X;’s,  i  =  M  +  I,  . . . ,  P  .It 
follows  from  (16)  and  (17)  that 

[Bgj  Pnj]  ^ (ifSc7)  [Bsj  ^Nj]  —  [^sj^Nj] 

X  ( A+P+ A  + )  [  Bgj  Pnj]  +  —  diag  [  A2] , 

or 

A^R^A^Bsj  =  Bs^Ai  -  B^jic^D 

=  Bg^^Diagl/Ui-  •  -MmI  (16) 

and 

Bf3,,A,B.A«  =  A2Bg,,-(cr^/)Bg,,  =  0.  (19) 

SinceP+A?  isofrankM,  it  follows  from  ( 19 )  that  the 
{(P  -  M)  X  M)  matrix,  BnjA+  must  satisfy 

B»^A+  =  0  or  ATb;^.;  =  0  (20) 

where  the  star  superscript  denotes  the  operation  of 
complex  conjugation.  Since  B^jBgj  —  I,  (18)  implies 

(Bs,^A+)P+(A  +  Bs,,)  =  Diagl^i  •  •  • 


which  in  turn  implies  that  the  (M  X  M)  matrix 
A  +  Bgj  is  nonsingular. 

It  then  follows  from  ( 18)  that 

=  Bs.,Diag[M,-  •  + 

“  ^  sjC+,  ( 21 ) 

where  =  Diag[Mr  ’  'Ma/] +  is,  obvi¬ 

ously,  nonsingular.  The  role  of  C+  in  ( 21 )  is  analogous 
to  the  role  ofCin  (14).  Since  =  A“,  (13)  leads  to 

AVBn  =  0  and  =  0.  (22) 

Using  (14),  (21),  and  (22),  an  algorithm  for  esti¬ 
mating  e^"‘’s,  t  =  1,  2,  . . . ,  M,  is  derived  next. 


After  associating  a  companion  matrix  D  of  order  P 
with  the  polynomial  on  the  left-hand  side  of  (27), 


0  1  0  •••  0  ~ 

0  0  1  •  ■  •  0 


D  = 


,  (28) 


0 

0 


0 

bo 


b 


P-1 


0 

b. 


b 


P-1 


1 

^P-2 

^P-1 


it  is  clear  that  the  polynomial  equation  in  ( 27 )  may  be 
rewritten  as 


det  [2/  —  D]  =  0. 


2.2.  Derivation  of  the  Algorithm  for  Estimating 
Parameters 
After  defining, 

aiz)  -  [1  2*  •  •  (23) 

it  is  clear  from  (22)  that  2^  -  t  =  1,  2, .  .  .  ,M,  are 
the  M  roots  of  the  ( P  —  1 )  th-degree  polynomial  equa¬ 
tions 

[a{z)]^b  =  60  +  6i2  +  622^ 

+  •  •  •  +  6p_i2^“^  =  0,  (24) 

where  P  =  [  60  ^1  *  *  *  bp„i  ]  is  any  one  of  the  columns  of 
Pn  or  B^j.  Therefore, 


The  desired  roots,  2,  =  f  =  1,  2,  .  . .  ,  M,  on  the 
unit  circle  are  also  the  M  unit  magnitude  eigenvalues 
Zi,i  —  1,2,.  .  .  ,  Mof P.  For  each  such  eigenvalue,  the 
associated  eigenvector  is 

a(2,)  -  [1  2,27-  •  -zrM'.  (29) 

Therefore,  the  following  eigenequations  hold. 

Da{Zi)  ^  Zia(Zi),  f  =  1,  2,  .  .  .  ,  M.  (30) 
The  set  of  equations  in  ( 30 )  may  be  rewritten  as 
P[a(2i)-  •  -0(2^)]  =  [a(2i)*  •  •a(2M)]A^,  (31) 

where 


Alb  =  0. 


(25) 


=  Diag[2i  ^2  -  ‘ 


Each  column  of  and  is  associated  with  an 
equation  similar  to  (24).  The  set  of  these  2(P  -  M) 
equations  can  be  written  as 


Since  2,-  -  ^  1,  2,  .  .  .  ,  M, 

-  [q(2i)  0(22)  •  *  *0(2^)], 


(32) 


(33) 


[ci{z)Y[bi-  •  'b^{p~M)]  —  0, 

and  from  (25) , 


^  +  [^1*  ’  ’^2(p-M)]  “9-  (26) 

Without  loss  of  generality,  it  may  be  assumed  that 
0  in  (24),  for  if  =  0,  the  monomial  of 
highest  degree  is  selected.  Multiply  (24)  by  2,  divide 
throughout  by  6p_i ,  and  rearrange  as  follows: 


Op_l 


+ 


b. 


2  =  0.  (27) 


and,  therefore,  (31)  is  expressible  as 


DA^  -  A^A„  (34) 

where  A^  is  the  diagonal  matrix  containing  the  M  so¬ 
lutions  desired.  The  theorem,  stated  and  proved  be¬ 
low,  shows  how  information  from  the  signal  and  noise 
subspace  eigenvectors  provides  an  estimate  of  the  pa¬ 
rameters.  Before  applying  the  theorem,  it  is  necessary 
to  construct  from  the  standard  model  of  the  ( P  X  p) 
correlation  matrix  S,  the  (P  X  Af )  eigenvector  matrix 
Bs.  whose  columns  span  the  signal  subspace,  and  the 
(P  X  (P  -  M))  eigenvector  matrix  whose  col¬ 
umns  span  the  noise  subspace.  Then  the  (P  X  P) 
companion  matrix  D  shown  in  (28)  is  constructed 


from  Bn.  while  the  (M  X  M)  matrix  (JBs)“B(JBs) 
is  generated  from  Bs  and  D . 

Theorem  1.  The  (M  XM)  matrix  ( JBs)^D(JBs) 
is  nonsingular,  whose  eigenvalues  uniquely  estimate 
the  M  distinct  parameters  e-'"*,  i  =  1,  2, . . . ,  M,  in  the 
standard  model  of  the  signal  correlation  matrix.  It  also 
holds  for  BsDB%. 

Proof.  The  matrix  D  in  (28)  can  be  written  as 


where  I  is  the  identity  matrix  of  order  (P  -  1 )  and  / 
is  a  row  vector  having  (P  —  1)  elements. 

A+  can  be  decomposed  into 


\h 


where 


and  h2=[l---l].  (35) 


Using  (21), 


B^jDBsj=  (C;')“a!?da+c+' 


0  I]\h2 


=  {Cl^)^[A^A^  +  hfg^A2]CZ^. 
From  (35)  it  follows  that  A2  =  where 
=  diag[e^“‘  •  •  •  e-'""]. 
Equations  (25),  (27),  and  (28)  imply  that 


1]A+  =  [-g^  1] 


0  or  g^Ai  =  hi- 


Therefore, 


BhDBsj=  + 

=  C+A,C;S  (37) 

since  (C;M«A?A,  =  =  C,. 

Thus,  we  see  that  B^jDB^j  =  C+AeC+^  holds  ex¬ 
actly,  and  since  C+  is  nonsingular,  the  eigenvalues  of 
BsjDBsj  uniquely  correspond  to  the  diagonal  ele¬ 
ments  in  A^  which  are  the  desired  parameters.  It  is 


straightforward  to  see  from  ( 14 )  and  A+  -  A  *  that  the 
{MX  M)  matrix  B^DBl  generated  by  replacing  B^j 
with  Bt  in  ( 37 )  is  also  nonsingular,  whose  eigenvalues 
are  i  =  1,  2,  . .  . ,  M. 

The  eigenvectors  of  the  (M  X  M)  matrix 
iJBs)^D{JBs)  which  can  be  used  in  the  estimation 
of  or  R  may  be  related  to  the  matrix  C+  in  ( 21 )  in  a 
manner  stated  and  proved  in  Theorem  2  below. 

Theorem  2.  The  (M  X  M)  matrix  E  of  the  M  ei¬ 
genvectors  of  the  matrix  (JB^)^D(JB^),  determined 
subject  to  the  constraint, 

(JBs)iE=[l  1--1].  (38) 

where  (JBq)i  denotes  the  first  row  of(JBQ),  is  the  ma¬ 
trix  C+  satisfying  the  relation,  A+  =  BsjC+  JBsC+ 
in  ( 21). 

Proof.  Let  2^ ,  i  =  1,  2,  . . .  ,  M,  represent  the  dis¬ 
tinct  eigenvalues  of  ( JBs )  ^ D  {JBg ) .  Since  ( JB^ )^D- 
(c/jBg)  is  nonsingular,  therefore,  the  matrix 

[{JBs)^D{JBs)  -  zj]  is  of  rank  M-  1,  and  nullity 

1.  Consequently,  the  eigenvector  c,  in  the  eigenequa- 
tion, 

[(JBs)»D(JBs)  -  2./]S,  =  0,  (39) 

has  only  one  nonzero  free  variable  u, .  Therefore,  for  i 
=  1,  2,  . . . ,  M 

e,  =  {e!)u,,  (40) 

where  ef  is  an  (M  X  1)  vector  of  constants.  Conse- 
quently, 

E  =  [e^  €2*  •  * 

=  e|  •  •  •  eM]diag(UiU2-  •  •  u^)-  (41) 

Since  E  is  constrained  as  in  (38) ,  therefore  following 
the  substitution  of  (41)  into  (38)  it  follows,  for  i  =  1, 

2,  . . . ,  M  that 


'  WBshef  • 

Since  the  ufs  are  unique,  E  is  also  uniquely  deter¬ 
mined.  If  (JBs)i  denotes  the  ith  row  of  JBg  for  i  =  1, 
2, ...  ,P  then  from  ( 21 ) 


{JBs)rC^ 

{JBs)pC, 


Note  that  each  element  in  the  first  row  of  A+  must  be 
1.  To  satisfy  this  condition,  replace  C+  by  £1  in  (43) 
and  apply  the  constraint  in  (38)  to  get 


A, 


1  1 

WBshE 

{JBs)pE 


(44) 


Clearly,  then,  E  is  the  unique  matrix  satisfy¬ 
ing  (21). 


Equation  (26)  provides  2(P  —  M)  eigenequations 
through  Theorem  1  to  decide  only  one  set  of  eigen¬ 
values.  There  may  be  several  ways  for  solving  the 
overdetermined  eigenvalue  problem.  In  the  example 
given  below,  simply  the  mean  of  the  2  (P  -  M)  coeffi¬ 
cient  sets  in  (26)  was  used  as  b  to  give  one  eigenequa- 
tion. 


which  for  N  consecutive  received  samples  leads,  im¬ 
mediately,  to 

Xi(n)  ^  [xi(n)  Xi(n  +  !)•  •  ■  x^in  +  N  -  1)]‘ 

=  As(n)  +  v(n),  (46) 

where 


s(n)  =  [e  ^'“'^(n)  e  ■'“^^(n)-  •  • 
v{n)  =  [u(n)  v{n+ 1)  u{n  + N— !)]'■. 


A  = 


1  ■  •  •  1 


g  -y(iv-l)w;v 


Therefore,  subject  to  the  assumptions  made 


^^USTRATIVE  EXAMPLE  AND  COMPARISONS 

The  example  discussed  here  is  the  one  considered 
in  [6],  where  the  MUSIC  algorithm  was  applied  to 
estimate  the  parameters  of  a  real  sinusoidal  signal 
embedded  in  zero-mean  white  Gaussian  noise  under 
various  signal-to-noise  ratios  ( SNRs) .  It  becomes  es¬ 
sential  to  show  how  this  type  of  problem  may  be  fitted 
to  the  model  in  ( 5 )  so  that  the  technical  devices  of  the 
previous  section  may  be  applied.  Let  the  number  of 
complex  sinusoids  be  M.  Then,  the  received  signal  at 
one  sensor  is 


M 

^i(n)  =  2  +  v^{n),  (45) 

k^\ 

where  each  discrete  sinusoid  whose  angular  frequency 
is  required  to  be  estimated  is  represented  by 
The  component  s^(n),  comprised  of  the 
magnitude  and  phase  random  variables,  is  stochastic. 
The  magnitude  and  phase  are  characterized  by  proba¬ 
bility  distribution  functions  which  permit  the  process 
to  be  viewed  as  zero-mean  wide-sense  sta¬ 
tionary.  In  (45),  Viin)  is  a  random  variable  which  is 
associated  with  a  zero-mean  wide-sense  stationary 
uncorrelated  random  process  having  a  variance  . 
Equation  (45)  may  be  recast  in  the  form 

Xi(n)  = 


X 


Siin) 

s^in) 


+  Viin), 


Sm(^) 


E[x^{n)xf{n)]  =  E[iAs{n)  v{n)) 

X  (As(n)  +  y(n))^] 

=  AEls{n)s^(n)]A^^a^L  (47) 

The  correlation  matrix,  E[xi{n)xf{n)],  on  the  left- 
hand  side  of  (47)  may  be  estimated  by  the  M  X  M 
matrix 


5  =  (48) 

where  for  the  N  samples  of  the  received  signal, 
{xi(0)  Xi(l)  -  •  -Xi(iV-l)},the(iV-P+l)XP 
matrix  X  is  given  by 

■  x(P-l)  x(P-2)  •••  x(0) 

x(P)  x(P-l)  •••  x(l) 

x{N-l)  x{N-2)  x{N-P) 

To  the  problem  solved  in  [6],  we  apply  here  not 
only  the  new  method  proposed  through  Theorem  1 
but  also  another  subspace-based  method  called 

GEESE.  It  is  also  pointed  out  that  the  estimates  of 
the  parameters  in  [  6  ]  are  inferred  from  a  graph.  Here, 
we  derived  the  estimates  by  applying  the  Root- 
MUSIC  algorithm  [  7  ]  in  addition  to  the  GEESE  algo¬ 
rithm.  The  Root-MUSIC  method  requires  the  evalua¬ 
tion  of  the  roots  of  the  equation 

G  (2;)Pj^P^a(2:  ^)  =  0,  (50) 


where 


aiz)  =  [1  z*  •  • 


(51) 


The  correlation  matrix  S  is 


The  GEESE  method  [1,  pp.  47-76]  requires  the  com¬ 
putation  of  the  singular  values  of  the  matrix  pencil 
{Bi,  B2},  with 

[Bgi  •  *  *  ^s(p-i)  (52a) 

B2—  [Bq2  ^S3*‘*^sp]*>  (52b) 

where  Bqi  denotes  the  ith  row,  i  =  1,  2, . . .  ,  P,  of  the 
(P  X  M)  signal  subspace  eigenvector  matrix  Bq. 

The  problem  considered  in  [6]  involves  the  re¬ 
ceived  data  %  ( n ) ,  n  =  0,  1,  . . . ,  6,  modeled  by 

x{n)  =  +  + 

1/2  i2 

The  parameters  ooi  and  0^2  are  to  be  estimated  from  the 
data  matrix 

'  x{3)  x{2)  3c(l)  x(0)  ■ 

x(4)  x(3)  x(2)  x(l)  / 2  V 

x(5)  x(4)  x(3)  x(2)  • 

x(6)  x(5)  x(4)  x(3) _ 

The  correlation  matrix  S  is  generated  from  X  through 
the  equation, 

S  =  X”X. 

Three  cases  are  considered.  The  data  matrices  are 
given  for  SNR  values  of  10,  20,  and  30  dB.  The  esti¬ 
mates  for  wi  and  0)2  are  obtained  by  applying  the  new 
algorithm  presented  here,  followed  by  Root-MUSIC 
and  GEESE  algorithms.  The  correct  values  of  the  pa¬ 
rameters  chosen  in  the  simulation  were 

coj  =  —1  and  «2  =  -tl.  (54) 

Consequently,  roots  Zi  and  Za  in  the  z-plane  are  on  the 
unit  circle,  |z|  =  1,  where  Z;  =  exp  (jco,- ),  i  =  1, 2.  Since 
the  data  is  noisy,  the  estimated  roots  need  not  lie  on  a 
circle  of  unit  radius.  Therefore,  the  radii,  r j  and  r 2 , 
each  of  which  ideally  should  be  unity,  as  well  as  the 
angular  frequencies  are  given  for  each  of  the  three 
methods  applied  to  solve  the  problem. 

Case  1  ( SNR  =  10  dB ) .  The  data  matrix  X  is 

-0.99313  0.42526  1.82929  1.31736 

-0.90984  -0.99313  0.42526  1.82929 

-0.87901  -0.90984  -0.99313  0.42526 

-0.21668  -0.87901  -0.90984  -0.99313. 


2.63373  1.47147  -1.13353  -3.13129 

1.47147  2.76762  2.05893  -0.77045 

-1.13353  2.05893  5.34126  3.66901 

-3.13129  -0.77045  3.66901  6.24889. 


The  eigenvector  matrix  associated  with  S  is 

-0.47835  0.68927  0.40616  -0.36210 

0.71288  0.09170  0.69515  0.01253 

-0.49272  -0.35650  0.54186  0.58010 

0.14212  0.62402  -0.24122  0.72953. 


The  eigenvalues  of  S  are 

0.20357  0.58090  5.49965  10.70738. 

The  estimates  of  a)i  and  CO2  as  well  as  the  values  of  Tj 
and  r2,  each  of  which,  ideally,  should  be  unity  are 
given  below  for  each  of  the  three  algorithms. 


Algorithm 

ri 

^2 

^2 

New 

Root-MUSIC 

GEESE 

0.90905 

0.85560 

0.87124 

-1.00776 

-0.90655 

-0.90967 

0.90905 

0.85560 

0.87124 

1.00776 

0.90655 

0.90967 

Case  2  ( SNR  =  20  dB ) .  The  data  matrix  X  is 

-0.94613  0.40878  1.50696  1.14561 

-1.24504  -0.94613  0.40878  1.50696 

-0.68038  -1.24504  -0.94613  0.40878 

0.45395  -0.68038  -1.24504  -0.94613. 

The  correlation  matrix  S  is 

3.11427  1.32945  -1.85619  -3.66774 

1.32945  3.07530  2.25433  -0.82270 

-1.85619  2.25433  4.88332  3.13361 

-3.66774  -0.82270  3.13361  4.64561. 

The  eigenvector  matrix  associated  with  S  is 

-0.50255  0.62403  0.34020  -0.49224 

0.65479  0.13263  0.74396  0.01380 

-0.55649  -0.28544  0.53005  0.57261 

0.09496  0.71521  -0.22324  0.65546. 


The  eigenvalues  of  S  are 


0.01975  0.04224  5.53627  10.12025. 


The  estimates  of  coj  and  0)2  as  well  as  the  values  of  Tj 
and  rj  are  given  below  for  each  of  the  three  algo¬ 
rithms. 


Algorithm 

0)1 

^2 

0)2 

New 

0.97229 

-1.00502 

0.97229 

1.00502 

Root-MUSIC 

0.94576 

-0.97266 

0.94576 

0.97266 

GEESE 

0.94879 

-0.97364 

0.94879 

0.97364 

Case  3  ( SNR  =  30  dB ) .  The  data  matrix  X  is 

-0.93126  0.40357  1.40503  1.09130 

-1.35104  -0.93126  0.40359  1.40503 

-0.61757  -1.35104  -0.93126  0.40357 

0.66602  -0.61757  -1.35104  -0.93126. 

The  correlation  matrix  S  is 

3.51753  1.30539  -2.17842  -3.78401 

1.30539  3.23682  2.28371  -0.83815 

-2.17842  2.28371  4.82955  2.98271 

-3.78401  -0.83815  2.98271  4.19516. 

The  eigenvector  matrix  associated  with  S  is 

-0.50987  0.58999  0.32172  -0.53707 

0.63320  0.14285  0.76061  0.01142 

-0.57720  -0.26445  0.52162  0.56992 

0.07703  0.74938  -0.21420  0.62178. 

The  eigenvalues  of  S  are 

0.00196  0.00369  5.59116  10.18224. 

The  estimates  of  coi  and  <j)2  as  well  as  the  values  of  Cj 
and  r2  are  given  below  for  each  of  the  three  algo¬ 
rithms. 


Algorithm 

0)i 

^2 

CU2 

New 

Root-MUSIC 

GEESE 

0.99178 

0.98199 

0.98288 

-1.00201 

-0.99149 

-0.99182 

0.99178 

0.98199 

0.98288 

1.00201 

0.99149 

0.99182 

4.  CONCLUSIONS 


The  MUSIC  and  Root-MUSIC-based  algorithms 
explicitly  use  the  eigenvectors  spanning  the  noise 
subspace  to  estimate  the  direction-of-arrivals  and  the 
angular  frequencies  of  sinusoids.  Other  subspace 
methods  such  as  GEESE  and  ESPRIT  use  the  signal 
subspace  eigenvectors  in  the  parameter  estimation 
problems  mentioned  above.  The  present  method  uses 
both  the  noise  and  signal  subspace  eigenvectors  for 


the  stated  purpose  as  seen  from  Theorem  1.  Specifi¬ 
cally,  the  matrix  D  is  computed  from  the  noise  sub¬ 
space  eigenvectors  while  the  matrix  ( JRs)^D(  JBg) , 
whose  eigenvalues  have  to  be  computed,  also  depends 
on  the  signal  subspace  eigenvectors  explicitly. 

The  new  method  is  seen  to  give  better  estimates  for 
ct)i  and  0^2  than  Root-MUSIC  and  GEESE  for  the 
benchmark  problem  discussed  in  this  paper.  It  is  rea¬ 
sonable  to  expect  that  the  statistics  of  the  data  can 
influence  the  accuracy  of  the  results  derived  from  im¬ 
plementation  of  the  various  algorithms,  especially 
when  the  SNR  is  low.  However,  the  presented  method 
contributes  to  the  popular  subspace-theory-based 
techniques  through  the  development  of  a  new  algo¬ 
rithm  which  works  better  in  some  situations.  The 
conditions  on  the  data  under  which  such  improved 
performance  may  be  expected  and  optimal  solution 
schemes  for  the  overdetermined  eigenvalue  problem 
are  currently  under  study.  The  scopes  for  generalizing 
the  approach  to  the  parameter  estimation  problems 
for  two-dimensional  sinusoids  in  noise  is  also  worth 
investigating. 
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In  this  paper,  three  new  methods  called  the  PESS  meth¬ 
ods  are  presented.  The  first  is  a  1-D  based  2-D  PESS 
method  for  estimating  2-D  wavenumbers.  This  provides  the 
basic  framework  for  all  the  PESS  methods.  The  second  is 
the  Pairing-PESS  method  which  offers  efficient  pairing  in 
all  1-D  based  2-D  methods.  The  third  is  the  Direct  2-D 
PESS  method  which,  though,  in  the  category  of  direct  2-D 
methods  in  the  sense  that  it  is  not  based  on  any  1-D  signal 
processing  technique,  is,  nevertheless,  computationally  ef¬ 
ficient  even  when  compared  with  the  previous  1-D  based 
2-D  algorithms.  This  algorithm  is  different  from  the  other 
direct  2-D  algorithm  in  that  it  does  not  depend  on  spectral 
search  in  2-D  parameter  space.  Examples  are  presented  to 
illustrate  the  methods  and  a  comparison  is  made  with  the 
matrix  enhancement  matrix  pencil  method  recently 
proposed.  ©  1995  Academic  Press,  Inc. 


1.  INTRODUCTION 

Extensions  of  l-D  subspace  methods  with  the  ob¬ 
jective  of  estimating  the  wavenumber  of  a  2-D  sinusoi¬ 
dal  signal  in  noise  as  well  as  the  direction-of-arrival 
(DO A),  have  been  studied.  The  obstacle  to  finding 
computationally  efficient  high-resolution  algorithms 
in  such  cases  is  the  factorability  problem  associated 
with  bivariate  polynomials  in  ROOT-MUSIC  and 
ROOT-MIN-NORM-based  methods  (the  zeros  of  bi- 

*  This  research  was  supported  by  the  Office  of  Naval  Research 
under  Grant  N00014-92-J-1755. 
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variate  polynomials  trace  continuous  algebraic  curves 
[1];  see  [2]  for  the  testing  for  zeros  on  the  distin¬ 
guished  boundary  of  an  unit  polydisc  of  a  multivariate 
polynomial)  and  the  so-called  problem  of  pairing  bi¬ 
variate  polynomial  zeros  in  2-D  generalizations  of 
ESPRIT.  The  2-D  parameter  estimation  problem  re¬ 
ferred  to  above  has  been  tackled  by  two  approaches. 
The  first  approach  {1-D  based  2-D  approach)  requires 
the  decomposition  of  a  2-D  problem  into  two  separate 

1- D  problems  followed  by  the  construction  of  a  set  of 

2- D  solutions  from  the  set  of  1-D  solutions  by  pairing 
algorithms.  The  second  method  (direct  2-D  approach) 
solves  the  2-D  problem  directly  by  searching  for  spec¬ 
tral  peaks  or  manipulating  optimization  algorithms  in 
multidimensional  parameter  space.  Although  the 
former  is  attractive  computationally,  it  results  in  sub- 
optimal  solutions  and  suffers  from  the  difficulty  in 
implementing  pairing.  The  latter  can  be  modified  to 
produce  optimal  solutions  but  the  computation  cost  is 
very  high. 

A  computationally  efficient  1-D  based  2-D  wave- 
number  estimation  method  was  proposed  in  [5].  This 
method  identifies  the  state  equation  form  inherent  in 
the  2-D  covariance  data  matrix.  The  given  sample 
data  matrix  is  singular  value  decomposition-decom¬ 
posed  (SVD-decomposed)  to  yield  the  factors  of  the 
model.  This  model  gives  the  two  state  transition  ma¬ 
trices  in  diagonal  form.  The  diagonal  elements  of  each 
are  the  desired  1-D  signal  zeros  in  each  of  the  two 
dimensions.  This  method,  called  the  state-space 
method,  assumes  distinct  frequencies  in  either  di¬ 
mension  and  ignores  noise.  The  subsequent  pairing 
procedure  uses  the  amplitudes  of  the  signals.  A 
method  similar  to  the  state-space  method,  called  the 
matrix  approximation  method,  was  proposed  in  [6], 
where  the  original  covariance  data  matrix  is  recon¬ 
structed  in  a  least  squares  sense. 
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Separable  extensions  of  ESPRIT  to  the  2-D  case 
have  been  discussed  as  1-D  based  2-D  DOA  estima¬ 
tion  methods.  Procrustes  rotations  (PRO)-based  ES¬ 
PRIT,  proposed  for  compensating  imperfect  array 
data,  was  applied  separately  and  the  estimated  azi¬ 
muth  and  elevation  angles  were  paired  by  evaluating 
the  array  response  vectors  corresponding  to  ail  the 
possible  combinations  of  the  estimates  and  then 
searching  for  those  which  are  the  closest  to  the  signal 
subspace  [7].  In  [16],  an  implementation  of  ESPRIT 
was  given  based  on  delta  array  geometry,  where  phase 
differences  are  estimated  and  used  as  the  criteria  for 
pairing  the  parameters  computed  separately.  Assum¬ 
ing  an  array  of  sensor  triplets  with  arbitrary  displace¬ 
ments,  two  matrix  pencil  equations  were  set  up  and 
two  sets  of  1-D  DO  As  were  computed  in  [15]  using  the 
SVD  of  three  data  matrices.  Observing  that  pairing  is 
automatic  in  the  noise-free  case,  noise  cancelation 
was  tried  through  the  estimation  of  the  perturbation 
matrices  for  the  noisy  case.  An  extension  of  ESPRIT 
to  2-D  case  using  the  concept  of  marked  subspace  was 
proposed  in  [10].  This  algorithm  introduced  a  mark 
operator  which  assigns  different  amplitudes  to  the 
roots  in  each  1-D  solution  set  but  the  same  amplitude 
to  each  correct  pair  of  roots.  The  paired  solutions  are 
identified  by  finding  two  roots  associated  with  the 
same  amplitude. 

Prony’s  methods  for  estimating  2-D  sinusoidal  fre¬ 
quencies,  amplitudes,  and  phases  were  proposed 
[4,17].  Particularly  in  [4],  the  signal  model  is  identi¬ 
fied  to  be  expressible  as  a  recognizable  rational  func¬ 
tion  [1,  p.  269].  The  coefficient  fitting  is  performed 
from  the  best  reduced  rank  least-squares  approximate 
obtained  via  SVD.  In  [8]  a  2-D  wavenumber  estima¬ 
tion  technique  based  on  a  matrix  pencil  (MEMP)  was 
proposed.  The  pairing  procedure  is  similar  to  that 
in  [7]. 

The  Fourier  method,  the  maximum-likelihood 
method,  the  maximum-entropy  method,  the  autore¬ 
gressive  method  and  the  MUSIC  method  for  estimat¬ 
ing  wavenumbers  or  DOA,  based  on  multidimen¬ 
sional  search,  were  presented  in  [11-13],  2-D  Toeplitz 
Approximation  Method  (TAM)  followed  by  2-D  MU¬ 
SIC  was  studied  in  [14]. 

An  extension  of  subspace-fitting  method  to  the  2-D 
case  was  investigated  in  [9],  where  optimal  and  subop- 
timal  procedures  for  estimating  2-D  DOA  were  for¬ 
mulated.  It  was  shown  that  for  uncorrelated  and 
correlated  sources  these  methods  were  significantly 
better  than  MUSIC.  The  computation,  however,  can 
be  prohibitively  expensive  without  good  initial  guess 
of  the  estimates. 

In  this  paper,  three  new  methods  called  the  PESS 
methods  are  presented.  The  first  is  a  1-D  based  2-D 
PESS  method  for  estimating  2-D  wavenumbers.  This 


provides  the  basic  framework  for  all  the  PESS  meth¬ 
ods.  The  second  is  the  Pairing-PESS  method  which 
offers  efficient  pairing  in  all  1-D  based  2-D  methods. 
The  third  is  the  Direct  2-D  PESS  method  which, 
though,  in  the  category  of  direct  2-D  methods  in  the 
sense  that  it  is  not  based  on  any  1-D  signal  processing 
technique,  is,  nevertheless,  computationally  efficient 
even  when  compared  with  the  previous  1-D  based  2-D 
algorithms.  This  algorithm  is  different  from  the  other 
direct  2-D  algorithms  in  that  it  does  not  depend  on 
spectral  search  in  2-D  parameter  space. 

2.  PROBLEM  FORMULATION 

It  is  assumed  that  the  2-D  array  of  discrete-valued 
data,  x(ni,  can  be  modeled  as 

M 

x(ni,  ^2)  =  2  ^2)-  (2.1) 

The  noise  process  v{n^y  n^)  is  assumed  to  be  zero- 
mean  and  white  over  and  with  variance  <t^.  The 
coefficient  a,*  represents  the  amplitude  and  phase  of 
the  ith  2-D  sinusoid.  The  number  M  of  2-D  sinusoids 
may  be  estimated  by  thresholding  a  set  of  eigenvalues 
of  a  covariance  matrix  as  mentioned  later.  In  this 
paper,  it  is  assumed  that  M  is  known. 

Assume  that  the  two-tuple,  (cjij,  W2i),  which  repre¬ 
sents  the  i-th  wavenumber,  is  distinct  for  each  f  in  1  < 
i  <  M,  that  is, 

\^2i  ^2kl  ’  1^21  =  <^2kl  ’ 

or  ,  whenever  i  ¥=  k, 

^  ^2k) 

Note  that  either  ajj*  or  0^2*  tie  zero  for  any  L  With¬ 
out  loss  of  generality,  the  trivial  case  when  both  are 
zero  is  excluded.  Also  the  alias-free  constraint  is  im¬ 
posed,  that  is,  is  unique  over  {— tt,  tt]  for  m  =  1,  2 
and  1  <  M, 

The  problem  to  be  solved  is  formulated  as  follows: 
Given  the  Ni  X  N2  data  matrix 

’  x(0,  0)  x(0,  1) 

^  ^  x(l,  0)  x(l,  1) 

MN,  -  1,  0)  x{N^  -  1,  1) 

x(0,  N2-I)  ' 

•  •  •  x(l,  N2-I)  ^  (2.2) 

•••  x(Ni- 1,  N2- 1). 


estimate  the  M2-D  wavenumbers,  (wi,,  for  1  <  i  < 

M,  in  paired  form. 


x(ni,  nj) 


3.  DERIVATION  OF  SUBSPACES  FOR  2-D 
SINUSOIDAL  SIGNALS 


+  o(ni,  nj).  (3.1) 


3.1.  Correlation  Matrix 


The  measurement  model  of  Eq.  (2.1)  can  be  ex¬ 
pressed,  using  the  vector  notation,  as 


For  selected  integers  Pj  and  P2  such  that  M  <  Pi  <  Ni 
and  Af  <  Pj  <  N2,  the  data  x(ni  -t-  nj  -f-  ^),  for  0  <  Aij 
<  Pj  —  1  and  0  <  ^  P2  —  1,  may  be  stacked  as  a 
(Pi -Pa)  X  1  column  vector, 


[x(ni,  nj) 

[rii,  712  +  1 

x(ni,  112  +  P2  -  1)  +  1,  ^2  +  Pz  -  1) 


x(ni  -I-  1, 1X2) 
x(ni  -h  1, 112  -t-  1) 


xirti  +  Pi  -  1,  nj) 

x(iii  -I-  Pi  -  1, 112  -I-  1) 


x(ni  +  Pi  -  1,  nj  +  P2  -  1)]‘, 


(3.2) 


where  the  superscript  t  denotes  transpose.  Substituting  the  measurement  model  in  Eq.  (3.1)  for  each  element  of 
x(iii,  112),  described  in  Eq.  (3.2),  we  get 


x(ni,  n^)  = 


I  (rti+ 1 ) 


g  J<"1  A/<  « 1  '*2) 


j(wil(rti+Pi-’l)+«2in2) 


(<j  H  n  1 +a>2 1  ( /I2+ 1 )  ) 


g  J  « 1 +7=’i  - 1 ) +tU2A/'*2) 

1  +‘*>2Af<”2'^  1 )) 


gy(wil(ni+Pi— l)-fc*»2i(n2+l)) 


^y{a»iirti+t.>2i(?»2+p2“l)) 

1  (^*1  1 1  (^■*■^2“  1 )  ) 


('il  1) +<‘*2Jlrf<'»2+^2“  1)  ) 


I  ("l+Pj  “  1) +<*>2 1  (^'*"'^2”  O  ) 


+Pl“~  1)  +W2JirfXrt2'*'^2“"  1)) 


Oi 

<h 


[v(ni,  1x2) 
v(n.i,  n2  +  l) 


v{ni  +  1, 1x2) 

11(1X1  ■+•  1, 1X2  -I-  1) 


u(ixi  Pi  -  1, 112) 

0(1X1  -I-  Pi  —  1, 1X2  -t-  1) 


(3.3) 


0(1x1,  Hz  +  P2  -  1)  0(1x1  +  1, 1x2  +  P2  -  1) 


0(1x1  -h  Pi  -  1, 1x2  +  P2  -  1)]‘. 


Denoting  the  rightmost  vector  in  Eq.  (3.3)  by  v(ixi,  1x2),  the  equation  may  be  rewritten  as 


x(ixi,  1X2) 


1 

1 

•  •  • 

gj"uwr 

•  •  • 

eJ<-2x 

.  .  . 

gJ«2A/ 

gy*‘>2A/gJ"LA/ 

g/«2li^2‘”0 

•  •  • 

gy‘*'2A/<^2‘"l) 

^‘u>a  i(P2- 1  )g7«i  1 

gy‘‘»2A#<^2”l)gy‘*»iA/ 

•  •  • 

,y(«12ni  +<*>22'*2) 


02^- 


+  v(ixi,  1X2). 


(3.4) 


After  introducing  notations, 


A  =  [a(e-'"”,  e-'"’®)  •  •  •  aCe'"'", 

(3.6) 

Eq.  (3.4)  is  simply  written  as 

x(ni,  rtj)  =  A-s(ni,  n^)  +  v(ni,  n^),  (3.7) 

where  s{n-^,  n^)  is  easily  identified  by  comparing  Eq. 
(3.4)  and  (3.7). 

We  see  that  a(e-'"*‘,  e-'"^*)  in  Eq.  (3.5)  may  be  ex¬ 
pressed  as, 

a(e-^"'',  e-'"»‘)  =  a2(e-'“’‘)  ©  a,(e-'"“),  (3.8) 

where  0  denotes  the  Kronecker  product  (let  A  =  B® 
C,  where  B  is  of  order  p  X  q  and  C  is  of  order  s  X  ir, 
then,  the  order  of  A  is  (ps)  X  (qt)  and 
(^)(>(i-i)+m).(to-i)+n)  “  and 

ai(e-'"‘')  =  [1  e-'"”-  •  . 

a^(gy<^)  =  [1  gM/.  .  ,g>«5i<a,-i)]t  (3  9^ 

Denoting  the  (i,  /2)-th  element  of  the  M  X  M  matrix 
s(ni,  ra2)s"(ni,  n^)  by  So,(ni,  n^),  we  have 

^2)  =  (3.10) 

If  #  (jo^,  for  m  =  1  or  2,  we  see  that 

lim  2 

Nm-^co  n„-0 

<  lim  2 

^  lim  - - - - — 

N„-.oo  1  - 

,.  1  2 

<  hm  - - - - ^  =  0. 


If  =  u)^,  for  m  =  1  or  2,  we  see  that 

1  ,  1 

lim  —  2  =  lim  -5^  2  1  =  1. 

Nn-^co  fi„-0  N„^oo  m  n„-0 

Since  (oj^,  Wg*)  ^  (wi*,  0)2*)  whenever  i  ¥=  k,  it  follows 
that 

2^  iVi-i  N2-1 

lim  lim  2  2  =  ^n-^k)* 

N2-^co  iVj-fcoo  r»i-0  n^—O 

where  5  is  the  Kronecker  delta  function.  So,  we  may 
view  the  vector  s(ni,  ^2)  to  consist  of,  as  elements, 
uncorrelated  random  variables.  Since  the  data  has  fi¬ 
nite  support,  we  assume  that  the  source  correlation 
matrix  S  defined  by 


S  =  E(s(ni,  n^ls^lni,  n^)]. 


(3.11) 


is  nonsingular,  where  the  superscript  H  denotes  com¬ 
plex-conjugate  transpose.  This  assumption  holds 
when  we  have  moderate  amounts  of  data,  and  the  2-D 
wavenumbers  are  not  too  close.  The  correlation  ma¬ 
trix  R  of  x(7Xi,  1X2)  has  the  structure 

R  =  £:[x(ni,  n2)x"(ni,  ^2)]  =  ASA^  +  (3.12) 

where  A  and  S  are  defined  in  Eqs.  (3.6)  and  (3.11), 
respectively. 

3.2.  Estimation  of  the  Correlation  Matrix  R 

Construct  a  P^X  matrix  Xin^,  n^)  expressed  by 
Xirii,  n^) 

x(n„  n^)  x(ni,  n,  4-  1) 

^  x(ni  +  1,  nj)  x{n^  4-  1,  n,  -t-  1) 

.x(ni  -b  Pi  -  1,  ra,)  x(ni  Pj  1,  -b  1) 


x(ni,  /i2  -b  P,  -  1) 

x(ni  4  1,  n,  4-  Pj  -  1) 

•  •  •  x(ni  4-  Pi  -  1,  rij  4-  P,  -  1) 

We  know  that,  from  Eq.  (3.2)  and  (3.13), 

x(ni,  n,)  =  vec  X(ni,  n^. 


(3.13) 


(3.14) 


where  the  vec  operator  stacks  the  columns  of  a  matrix 
into  a  vector.  (If  H  =  [hjh,-  •  -  hi,],  then  vec  H  = 
[hjhj  •  •  •  h^]‘.)  Actually,  x(ni,  n,)  is  a  window  through 
which  we  look  at  the  data  matrix  X  of  Eq.  (2.2).  The 
window  of  size  Pj  X  Pj  moves  horizontally  with  in¬ 
dexing  and  vertically  with  indexing  n,  within  the 


data  matrix  X.  Each  windowed  submatrix  X(ni,  of 
X  is  mapped  to  x(ni,  n^)  by  Eq.  (3.14).  In  fact,  vec 

X(nj,  nj),  for  0  <  ni  <  Nj  -  Pi  and  0  <  na  <  Nj  -  Pj,  is 
interpreted  to  be  the  realization  of  the  random  vector 
x(ni,  n^).  Define  the  reorganized  data  matrix  X,  to  be 

Xr  =  [vec  X(0, 0)  vec  X(l,  0) 

vec  X(0, 1)  vec  X(l,  1) 


vec  X(0,  N2  -  Pi)  vec  X(l,  -  P2) 

vec  X(Ni  -  Pi,  0) 

•  •  •  vec  X{Ni  -  Pi,  1) 


•  •  •  vec  X(iVi  -  Pi,  N2  -  P2)]  ,  (3-15) 

obtained  by  stacking  vec  X(rai,  n^),  for  0  <  n-i  <  Ni  — 
Pi  and  0  <  ^2  <  iVj  -  P2,  as  columns.  Clearly  X^  has 
P1P2  rows  and  [{N-^  -  Pi  +  1)(N2  -  P2  +  Dl  columns. 
Let  Jbe  the  permutation  matrix  of  appropriate  order 
defined  below. 

-0  •  •  •  0  1" 

0  •  •  •  1  0 

€/  —  ,  •  ... 

.i  •  6  6. 

By  analyzing  «/(x(ni,  n^))*  similarly  to  the  case  of 
x(ni,  n^),  discussed  in  the  previous  subsection,  we  can 
show  that  JX*  can  be  considered  to  be  a  realization  of 
x(ni,  n^).  Assuming  that  x(ni,  n2)  is  correlation  ergo- 
dic,  the  estimate  R  of  the  correlation  matrix  R  is  com¬ 
puted  from 


R  = 


2(lVi  -  Pi  +  l)(iV2  -  P2  +  1) 


(X,X"  + JX:XJJ). 

(3.16) 


From  Eqs.  (3.12)  and  (3.16),  since 


Since  the  first  Pi  rows  of  A  are  given  by 
[a,(e'"”)  •  •  •  ai(e''"'")],  and  Px>  M,  A  is  of  full  rank, 
provided  the  Wk’s,  for  1  <  i  <  M,  are  distinct.  When 
the  «i,’s  are  not  distinct,  consider  the  linear  equation 

M 

2  Cia(e-'“*%  e'"®*)  =  0 


which,  after  using  Eq.  (3.8),  becomes 

M 

2  Cja2(c-'"'‘)  ©  ai(e-'“’'‘)  =  0. 

t-1 

Suppose  that  all  the  column  vectors  of  A  in  Eq.  (3.6), 
whose  first  element  is  the  same  as  e^“'S  are  a(e^“'>, 
gj"j..),  for  1  <  <  m.  Then  the  above  equation  is  reor¬ 

ganized  as 


[2  Ci,a2(e-'"’“)]  ©  a,(e^'“")  -t-  2  e^""')  =  0, 


A- 1 


iGl?! 


where  wu  +  wu  for  any  i  £  S?i.  Since  ai(e-'"")  and 
ai(e'"’‘)  are  linearly  independent  for  i  £  7i,  and  ai(  • ) 
appears  as  in  Eq.  (3.17),  it  must  follow  that  Z^_i  c,^ 
a2(e'"*‘‘)  =  0,  which  in  turn  implies  that  c,^  =  0  for  1  < 

<  m,  due  to  the  fact  that  the  wji.’s  are  distinct  from  the 
assumption  that  (wi,-,  0)2,)  ^^2*)  whenever  i  k. 

Continuing  the  above  argument  for  each  e-'"'s  we  get  c, 
=  0  for  1  <  I  <  M.  Thus,  all  the  columns  of  A  are 
linearly  independent  and  so  A  is  of  full  rank,  if  (wi., 
wji)  (wu,  <*>2*)  whenever  i  #  k,  which  is  consistent 
with  the  results  of  [4]. 

Since  the  correlation  matrix  R  is  Hermitian,  and 
ASA^  has  rank  M,  the  eigendecomposition  of  R  has 
the  structure  described  by  [19,  p.  168] 

B“RB  =  B^ASA^B  + 

=  Diag[/ii  +  •  •  */*«+  a^‘  •  •  <t^],  (3.18) 


EIR]  =  R, 

the  estimator  in  Eq.  (3.16)  is  unbiased. 

3.3.  Derivation  of  Subspaces 

The  (P1P2)  X  M  matrix  A  in  Eq.  (3.6)  has  columns, 
a(e'"“,  e'"*0  -  a2(e''"*‘)  ©  ai(e'"’0.  for  1  <  i  <  M,  each  of 
which  can  be  expressed  by 


a(e>"i.,  e>"«)  = 


ai(e-'"'0 

e^"’‘ai(e-''">0 


where  B  is  a  unitary  matrix  of  order  P1P2  and  /i.  >  0, 
for  are  the  eigenvalues  of  ASA^.  The  num¬ 

ber  M  of  2-D  sinusoids  may  be  estimated  by  threshold¬ 
ing  the  sequence  of  eigenvalues.  Partition  the  unitary 
matrix  B  =  [bib2  •  •  ’bpiPj]  i^fo  B  =  [BgBj^],  where  Bs 
=  [bt  •  •  •  b*,]  is  the  matrix  of  eigenvectors  associated 
with  the  eigenvalues,  Mi  +  for  1  <  i  <  M,  and  B;v  = 
[bjif+i  •  •  •  bp,pj  is  the  matrix  of  eigenvectors  associated 
with  the  eigenvalues  (p-.  Through  the  same  argument 
as  in  [3],  Eq.  (3.18)  implies  that 


gMi(i’2-i)ai(e-'">0 


(3.17) 


A"B^  =  0, 


(3.19.a) 


and 


A  -  BsZ  (3.19.b) 

where  T  is  a  nonsingular  matrix.  These  results  will  be 
used  in  all  the  PESS  algorithms. 

4.  DERIVATION  OF  SUBSPACES  FOR 
SEPARATE  1-D  SIGNALS 


4.1.  Correlation  Matrices  and  Their  Estimation 

x(ni  +  ki,  rii),  for  0  <  Pi  —  1,  are  mapped  into  a 
Pi  X  1  column  vector 

Xi(ni,  n^) 

=  [xin^y  n^)x(n^  +  1,  ^2)  •  •  •  x{n^  +  Pi  -  1,  n^)Y,  (4.1) 

Substituting  the  measurement  model  of  Eq.  (3.1)  for 
each  element  of  Xi(ni,  722),  we  get 


X 


2 

me^i 


2 


+  Vi(ni,  n^).  (4.4) 


Define  Ai  =  [ai(e’^"'‘0  •  •  •  ai(e-^‘*''‘*«)].  Using  ai(  • )  given 
by  Eq.  (3.9),  Eq.  (4.4)  is  written  as 


Xi(ni,  ^2)  =  Ai  •  Si(ni,  02)  +  Vi(ni,  712),  (4.5) 

where  Si(ni,  ^2)  is  identified  by  comparing  Eq.  (4.4) 
with  Eq.  (4.5).  The  /zth  element  of  the  M^X  I  vector 
Si(ni,  n2)  can  be  rewritten  as  e-'"^**'*'* 

Also  the  elements  of  the  set  for  1  ^  <  Mi,  are 

distinct.  Thus,  using  similar  arguments  as  in  Sec,  3.1, 
we  can  conclude  that  asymptotically,  Si(ni,  n2)  is  a 
random  vector  with  a  nonsingular  correlation  matrix 

Si  =  P[Si(n„  n2)sf(ni,  n^)]. 


Xi(ni,  nj) 


’Ctl" 

_ 

gy(«ii(ni  +  l)+cJ2in2)  ,  ,  .  + 

(h 

^;(wn(ni+Pi-l)+W2ifl2)  .  ,  • 

+  [y(ni,  nj)  o(ni  +  1,  na)  •  •  •  o(ni  +  Pi  -  1,  raj)]'- 

(4.2) 

Denoting  the  last  term  of  Eq.  (4.2)  by  Vi(ni,  n^)  and 
reorganizing  the  first  term,  we  have 


1 

1  ■ 

= 

. 

+‘»>22"2) 

+  Vi(ni,n2).  (4.3) 


We  know  that  Wi^,  for  may  not  be  distinct. 

Suppose  that  for  1  <  <  Mi,  are  distinct  such 

that  the  two  sets,  {a)i,„}"_i  and  {wi,-^}*l.*i,  are  the  same 
after  ordering,  and  for  m  £  7*,  where  2^'i 

(number  of  elements  in  ?*)  =  M.  Note  that  one  of  the 
Will’s  may  be  zero.  Then  Eq.  (4.3)  can  be  rearranged 
into 


Xi(ni,  rtj) 


1 


1 


The  correlation  matrix  Pi  of  Xi(ni,  n^)  has  a  struc¬ 
ture  similar  to  that  of  P  in  Eq.  (3.12). 

Pi  =  P[Xi(n.i,  ?22)x^(ni,  JT^)]  —  AiSiA^^  4-  <P/.  (4.6) 

In  Eq.  (4.6),  Ai  and  Si  are  the  counterparts  of  A  and  S 
in  Eqs.  (3.6)  and  (3.11).  Actually,  Xi(ni,  7x2)  is  a  win¬ 
dow  through  which  we  look  at  the  data  matrix  X  of  Eq. 
(2.2).  The  window  of  size  Pi  X  1  moves  horizontally 
with  indexing  and  verticedly  with  indexing  tZi 
within  the  data  matrix  X  in  Eq.  (2.2).  Xi(ni,  712),  for  0 
<  rii  ^  Ni  —  Pi  and  0  <  712  <  N2,  are  interpreted  to  be 
the  realizations  of  the  random  process  Xi(7ii,  712).  De¬ 
fine  the  reorganized  data  matrix  X^^i  by 

[Xi(0,0)  xi(l,0) 

Xi(0, 1)  Xi(l,  1) 

X,i  = 


Xi(0,  N2  -  1)  Xi(l,  N2  -  1) 

Xi(iVi-Pi,0) 

Xi(iVi-Pi,l) 

•••  Xi(iVi-Pi,iV2-l)]  ,  (4.7) 

which  has  the  order  Pi  X  ((iVi  -  Pi  +  l)iV2).  Assuming 
correlation-ergodicity  of  Xi(72i,  7I2),  the  estimate  Pi  of 
the  correlation  matrix  Pj  is  computed  by 


Next,  the  data  elements x(ni,  02  +  fh),fotO  ^ 

—  1,  are  stacked  into  a  P2  X  1  column  vector 

=  [x(nx,  ^  +  1)  *  •  •  x(ni,  n,  +  P2  —  !)]'•  (4.9) 

Through  the  same  argument  as  in  the  case  of  Xi(ni, 
n,),  we  get 

X2(ni,  02)=  A^' 32(^1,  n^)  +  V2(ni,  rij),  (4.10) 

where  =  [a2(e-'"’‘')  •  •  •a2(e-'“®**«)],  M2  is  the  number 
of  distinct  W2/3,  and  defined  by  P[s2(nj, 
rii)  ]  is  nonsingular.  The  correlation  matrix  R2  of  X2(ni, 
n2)  has  the  structure  expressed  by 

R2  =  ElXiirii,  n^x^irii,  Uj)]  =  A2S2A2  +  (4.11) 

where  Aj  and  S2  are  defined  analogously  to  A^,  Si,  or 
A,  S.  x^(ni,  ru^  is  a  window  through  which  we  look  at 
the  data  matrix  X  of  Eq.  (2.2).  The  window  of  size  1  X 
P2  moves  horizontally  and  vertically  with  indexing  ^2 
and  Hi,  respectively.  Define  the  reorganized  data  ma¬ 
trix  X,.2  by 

X,2=  [X2(0,0)  x:2(0,l)) 

X2(l,0)  X2(l,l) 

•  •  •  •  •  • 

•  •  •  •  •  # 

X2(iVi-l,0)  X2(iVi-l,l) 

^2(0,  N2  -  P2) 

X2( 


X 


X2(l,iV2-P2) 


X2(iVi-l.N2-P2)],  (4.12) 

which  is  of  order  P2  X  (iVi(iV2  -  P2  +  D).  Assuming 
correlation-ergodicity  of  X2(ni,  02),  the  estimate  P2  of 
the  correlation  matrix  R2  is  computed  by 


where  Si  and  S2  are  unitary  matrices,  and  Mii,  for  1  ^  i 
^  Ml,  and  for  1  <  i  <  Mg,  are  positive.  The  num¬ 
bers,  Ml  and  Mg,  of  distinct  1-D  parameters  may  be 
estimated  by  thresholding  the  sequences  of  eigen¬ 
values.  We  know  that  1  <  Mi  <  M  and  1  <  Mg  <  M.  In 
PESS,  Ml  and  Mg  are  not  required.  Partition  the  uni¬ 
tary  matrices  Bi  =  [biibjg*  • -bipj  and  Bg  = 
[b2ib22  •  *  *  bgp^]  into  Bi  =  [Bi^Bi^]  and  Bg  =  [Bg5B2^] , 
respectively,  where  B^s  =  [bu*  •  ’bijv^j]  and  Bos  ~ 
[bgi*  •  -bgjvfJ  are  the  matrices  of  eigenvectors  asso¬ 
ciated  with  the  eigenvalues,  for  1  ^  i  ^  Mi, 

and  fjL2i  +  for  1  <  i  <  Mg,  respectively.  Bi^  = 
[bi(A/j+i)  •  •  •  bipj  and  Bgjy  =  V^2{M2+i)  *  *  *  bgp^]  ^re  the 
matrices  of  eigenvectors  associated  with  the  eigen¬ 
values  of  Ri  and  Bg,  respectively.  Using  the  argu¬ 
ment  in  [3],  Eqs.  (4.14.a)  and  (4.14.b)  imply  that 


(4.15.a) 

II 

(4.15.b) 

A^Bgjv  “ 

(4.16.a) 

and 

Ag  =  B25T2, 

(4.16.b) 

where  Ti  and  Tg  are  nonsingular  matrices.  These  re¬ 
sults  will  be  used  only  in  the  1-D  based  2-D  PESS. 

5.  DERIVATION  OF  THE  1-D  BASED  2-D  PESS 
METHOD 

5.1.  Interpretation  of  the  Noise  Subspaces 

Using  Eq.  (3.19.a),  we  can  easily  set  up  the  straight¬ 
forward  2-D  extension  of  the  MUSIC  algorithm  as 
follows.  Find  (a)i,  ciJ2)’s  where  the  peaks  of  the  function 
below  occur. 


P(wi,  O^g)  = 


(5.1) 


4.2.  Derivation  of  Subspaces 

Since  the  correlation  matrices  Bi  and  Bg  are  Her- 
mitian,  and  AjSjAf  and  AgSgA^  are  of  rank  Mj  and 
Mg,  respectively,  the  eigendecompositions  of  Bi  and 
Bg  have  the  structures  described  by 

BfBiBi  =  BfAiSiAfBi  +  a^I 

=  Diag[)Ltii  +  0*2  •  •  •  MiAfi  +  0^  (4.14.a) 

B^BgBg  =  B^AgSgAfBg  + 

-  Diag[/X2i  +  •  M2M2  +  •  •  •  (T^],  (4.14.b) 


This  algorithm  suffers  from  prohibitive  computa¬ 
tional  load,  and  the  2-D  generalization  of  the  ROOT- 
MUSIC  algorithm  is  not  routine  because  the  zeros  of 
bivariate  polynomials  trace  continuous  algebraic 
curves  [1],  and  the  straightforward  extension  of  the 
Vandermonde  system  to  the  2-D  situation  cannot  be 
derived.  In  this  section,  an  indirect  2-D  generalization 
of  the  Vandermonde  system  is  defined  and,  based  on 
it,  a  new  method  for  estimating  the  wavenumber  pa¬ 
rameters,  (wii,  cjg,-),  for  1  <  i  <  M,  is  developed. 

For  each  column  bi;  of  Bi;^  given  by  Eq.  (4.15. a), 
construct  polynomial  equations 


b^'a,(z,)  =  0, 


(5.2.a) 

and,  for  Mj  +  1  ^  i  Pi, 

Jbi,ai(2i)  -  0,  (5.2.b) 

where  ai(2i)  =  [1  Zj-  •  Eq.  (5.2.a)  and  (5.2.b), 

for  Ml  +  1  <  i  <  Ply  include  the  separate  1-D  signal 
zeros,  for  1  <  ^  <  Mi,  one  of  which  may  be  e*'  *®,  as 

implied  by  Eq.  (4.15.a).  Using  Eqs.  (5.2.a)  and  (5.2.b), 
construct  the  set  of  (Pi  —  Mi)  polynomial  equations 
each  of  degree  Pi,  as  shown  below. 

I  .  z*  ^  .  (b"  +  Jbu)ai(z,)  =  0, 

k^O 

for  Ml  +  1  <  i  <  Pi.  (5.3) 

Let  the  resulting  set  of  (Pi  -  MJ  polynomial  equa¬ 
tions  be  described  by 

2  =  0,  1  ^  i  <  Pi  -  Ml  +  1.  (5.4) 

k-O 

The  (Pi  —  Ml)  sets  of  coefficients,  for  1  ^  i  < 

Pi  —  Ml  +  1,  are  interpreted  as  the  statistical  data  set 
for  a  generic  set  of  coefficients,  denoted  by 
Without  loss  of  generality,  we  can  assume  that  tip^ 

0.  If  tif,  =  0,  for  m  <  A*  <  Pi  and  ti(,„-i)  ^  0,  then  we 
multiply  both  sides  of  Eq.  (5.4)  by  Manipulate 

Eqs.  (4.16.a)  and  (4.16.b)  in  a  similar  way  to  get  a 
corresponding  set  of  equations,  given  next,  in  the 
other  dimension. 


b^^a2(z2)=0,  (5.5.a) 

Jh2i^{z2)  =  0,  (5.5.b) 


and,  for  M2  +  1  <  i  ^  P2, 


5.2.  A  2-D  Vandermonde  System 

A  2-D  Vandermonde  system  is  defined,  based  on  a 
set  of  1-D  Vandermonde  system  given  below.  Eq.  (5.4) 
and  (5.7)  can  be  reformulated  as  1-D  Vandermonde 
systems  expressed  by 


Di  •  ai(zi)  —  ai(zi)  •  Zi  (5. 8. a) 

D2  •  a2(z2)  “  3^(22)  *  Z2  (5.8.b) 

where  the  companion  matrices,  for  m  =  1  and  2,  are 


“  0  1  0 

0  0  1 


0 

0 


D 


m 


6  6  6 
^m.O  _  ___  ^m.2 


i 

(5.8.C) 


Taking  the  Kronecker  products  of  the  left-hand  sides 
and  the  right-hand  sides  of  Eq.  (5.8.a)  and  (5.8.b),  we 
get 

{D2 ©  Di)  •  [a2(z2)  0 ai(zi)]  =  [a2(z2)  0  ai(zi)]  •  (ziZg), 

(5.9) 


which  will  be  called  a  2-D  Kronecker  product  Vander¬ 
monde  system.  From  the  property  of  the  Kronecker 
product,  we  know  that  is  the  eigenvalue  as¬ 

sociated  with  a2(e^"^0  ©  ai(e-^"'0,  which  is  an  eigenvec¬ 
tor  of  {D2  ©  Di),  for  For  arbitrary  complex 

numbers  gi  and  ^2>  Eqs.  (5.8.a)  and  (5.8.b)  can  be 
rewritten  as 


*  ai(zi)  —  ai(zi)  •  (giZi)  (5. 10. a) 

ig2P>2)  •  a2(-^2)  =  ^(^2)  •  (^2^2)-  (5.10.b) 


2  47^^^  •  4  =  22  •  (b|^  +  Jb2,)a2(z2)  -  0.  (5.6) 

k^O 

Finally,  we  have  a  generic  polynomial  equation  repre¬ 
senting  Eq,  (5.6),  as  given  below. 


From  a  property  of  the  Kronecker  sum  denoted  by  ©, 
(B  ©  C  —  B  ©  /c  +  /b  ©  C,  where  Iq  and  Iq  are  the 
identity  matrices  of  the  same  order  as  C  andB,  respec¬ 
tively)  it  follows  that 


P2 

2  ^2k^  “ 


*-0 


with  tzp^  0. 


(5.7) 


[(£202)  ©  i£iDi)]  •  [a2(z2)  ©  ai(zi)] 

=  [a2(z2)  ©ai(zi)]  •  (^iZi  4-^222).  (5.11) 


We  know  that  Eq.  (5.4)  and  Eq.  (5.7)  include  the  sepa¬ 
rate  1-D  signal  zeros,  for  I  k  <  Mj,  and 
for  1  r  <  M2,  respectively.  When  Mi  and  M2  are  not 
available,  M  replaces  those.  In  simulation,  averages  of 
the  coefficient  sets  were  substituted  in  the  generic 
polynomial  equations. 


Equation  (5,11)  implies  that  is  the  ei¬ 

genvalue  associated  with  a2(e''‘^)  ©  ai(e'"^0,  which  is 
an  eigenvector  of  [(^2^2)  ®  (^1^1)  ]>  for  1  <  i  <  M, 
where  either  gi  or  g2  must  be  non-zero.  Equation 
(5.11)  is  called  a  2-D  Kronecker  sum  Vandermonde 
system. 


Multiplying  both  sides  of  Eq.  (5.9)  by  a  complex 
number  and  adding  the  resulting  equation  to  Eq. 
(5.11),  we  have 

[g^iD^  ©  DO  +  {(g^DO  ©  (SiDO}]  •  M^2)  ®  aiUi)] 

=  [02(^2)  ©  ai(2i)]  •  ig3ZiZ2  +  §2^2  +  (5-12) 

which  can  be  expressed  by 

Dig)  •  a.{Zi,  Z2)  =  \(g)  •  a(2i,  22).  (5.13.a) 

where 

D(g)  =  ^3(^2  ®  ^1)  +  {(^2^)  ©  isM 

©  DO  +  g^{D^  ©  IpO  +  giilp,  ®  DO, 

(5.13.b) 

Mg)  =  gzZl^2  +  g2^2  +  gl^U  (5.13.C) 

and  g  =  \gi,  g2,  gzV  is  a  non-zero  vector  (/^  denotes  the 
identity  matrix  of  order  0-  Eq.  (5.13a)  is  called  a  2-D 
Vandermonde  system.  We  know  that  + 

-H  is  the  eigenvalue  associated  with  a(e^"*S 
for  1  <  i  ^  M,  which,  for  any  non-zero  complex 
vector  g,  is  an  eigenvector  of  D(g).  Similarly,  the  M 
eigenvalues  are  denoted  by  Xi(g),  for  1  <  i  <  M.  Then 
from  the  definition  of  the  matrix  A  and  Eq.  (5.13.b), 
we  have,  for  any  nonzero  g, 

D(g)-A  =  A-A(g),  (5.14) 

where  A(g)  =  Diag[Xi(g)  •  •  •  XAf(g)]  and  Xj(g)  = 

g^e^<^u . 

5.3.  Obtaining  the  Parameter  Vectors 

Using  Eqs.  (3.19.b)  and  (5.14),  the  2-D  parameter 
vectors,  (a;i„  ^2^),  for  1  <  i  ^  M,  are  obtained  from  the 
results  of  the  following  theorem,  which  shows  the  sig¬ 
nal-selectivity  of  signal  subspace.  The  theorem  is  for¬ 
mulated  and  proved  in  a  general  framework. 

Theorem  5.1.  Consider  matrices,  D  of  order  p  X  p, 
A  and  E  each  of  order  pX  q  such  that  p  >  q,  let  E  be  of 
full  rank,  SpanfAj  =  Span(E),  and  DA  =  A  A,  where  A 
is  a  (qXq)  diagonal  matrix.  Then,  the  following  eigen- 
decomposition  exists 

[{E^E)^'^E^DE]  =  TAT“S  (5.15) 

where  A  is  the  same  as  A  after  column  or  row  permuta¬ 
tion.  Furthermore,  for  each  eigenvalue  X  of  multiplicity 
m  in  A, 

Span  (A  J  =  Span(£Tx),  (5.16) 


where  A^  and  are  the  matrices  of  m  eigenvectors 
associated  with  X  in  A  and  T,  respectively. 

Proof.  After  substituting  A  =  ET,  where  T  should 
be  nonsingular  matrix,  for  A  in  DA  =  AA,  the  least- 
squares  solution  for  TAT"^  is  obtained  as  {E^E)'^- 
E^DE  =  TAT“\  This  implies  that  {E^E)^^e^DE  is 
diagonalizable.  Pre-multiplying  and  post-multiplying 
both  sides  of  Eq.  (5.15)  by  E  and  T,  respectively,  we 
get,  E{E^E)^^E^DET  =  ETA,  Since  Span(ET)  = 
Span(A),  therefore,  ET  =  A Wfor  some  nonsingular 
matrix  W.  Therefore,  EiE^EJ'^^E^DAW  =  AWA, 
which  reduces  to  E(E^E)~^E^AAW  =  AWA,  since 
DA  =  A  A.  From  the  fact  that  AA  W  E  Span(E)  and 
E(E^E)~^E^  is  the  projection  operator  to  Span(E),  it 
follows  that  A  A  W  =  A  WA.  Since  A  is  also  of  full  rank, 
therefore,  AW  —  WA,  or  A  =  WAW^^.  Since  A  and  A 
are  diagonal  matrices,  and  the  last  equation  repre¬ 
sents  an  eigendecomposition,  therefore,  A  =  A  within 
permutation. 

Suppose,  for  convenience  of  argument,  that  A  =  A, 
and  Xj,  for  1  ^  i  <  r{^q),  are  the  distinct  elements  in  A, 
and  X,-  is  of  multiplicity  m^.  Then,  it  follows  that  A 
=  WAW~^,  where  A  =  Diag[Diag[Xi  •  •  •  XJ 
•  •  -DiaglX,.*  •  -XJ]  and  W  =  Diag[WiW’2*  •  •  WO, 
where  the  Wfs  are  nonsingular  matrices.  Thus,  from 
ET  ^  AW,  it  follows  that  A  =  ETW^\  or  A,^  = 
ET^WJ\  for  1  <  i  <  r.  ffl 

The  matrices,  A,  B5,  and  D  given  by  Eq.  (3.19.b) 
and  (5.14)  satisfy  the  assumption  of  Theorem  5.1.  If 
we  have  an  eigendecomposition, 

{B^Bs)-^B^D{g)Bs  =  T(g)A(g)T(g)-\  (5.17) 

where  A(g)  has  the  eigenvalues  X£(g)  of  multiple  m^, 
for  1  <  i  ^  r  with  r  <  M,  and  Tj(g)  is  the  matrix  of  all 
the  eigenvectors  corresponding  to  Xj(g),  then  we  have 
that 

Span(Ax.(g))  =  Span(B5T.(g)),  for  1  <  i  <  r.  (5.18) 

If  mi  =  1,  then  Eq,  (5.18)  implies  that 

==  BsTi(g)  •  (5.19) 

where  Tj(g)  is  actually  a  column  vector,  is  actu¬ 
ally  e^"^),  for  some  k,  and  is  a  scalar.  When 

we  scale  BsTi(g)  so  that  its  first  element  is  1,  let  the 
second  element  be  and  the  {Pi  +  l)th  element  be 
g>u>2A.  Then  {coik,  <*^2k)  is  a  desired  parameter  vector.  If 
m,-  =  1  for  all  i,  then  Eq.  (5.18)  implies  that 

A^BsT(g)^T,  (5.20) 

where  F  is  a  nonsingular  diagonal  matrix.  When  the 
scaling  is  done  so  that  the  first  row  of  BsT{g)  is 


and 


[1  •  •  •  1],  the  M  parameter  vectors,  (wj,-,  W2i),  for  1  <  i 
<  Af,  can  be  easily  identified  from  Eq,  (5.20). 

If  Mi  >  1,  we  apply  the  Theorem  5.1  again.  The 
matrices,  Ax.(g),  SsT^(g),  and  D(g')  given  by  Eq.  (5.18) 
and  Eq.  (5.14)  satisfy  the  assumption  of  Theorem  5.1. 
Note  that  g'  is  different  from  g,  so  that  the  eigen¬ 
values  of 

{Ti{g)^B^BsTi{g)n{g)^^^  (5.21) 

are  distinct.  Denote  by  Ti{g,  g')  the  m,-  X  eigenvec¬ 
tor  matrix  of  Eq,  (5,21),  then,  we  have 

Ax.,,)  =  BsT,(g,g^)-r,,  (5.22) 

where  is  a  non-singular  matrix.  Thus,  is  de¬ 

termined  uniquely  by  selecting  Fj  so  that  the  first  row 
ofEq,  (5.22)  is  [1-  •  •  1]. 

When  the  estimate  A  of  A  has  been  obtained,  each 
column  of  A  represents  a  paired  parameter  vector.  In 
order  to  compute  the  parameters,  the  structure  of  A  is 
exploited.  We  know  that 


1 

1  • 

1*4>2 

1*4>2*^i 

1  • 


(5.23) 


where  #j  =  Diag[^""  •  •  •  e'"*"]  and  $2  Diag 
Let  A  be  defined  by  the  column  vec¬ 
tors  named  below. 


•  •  •  (5.24) 

Then  it  follows  that 


Sia^i  =  and  =  S2b7  (5.25) 


Let  the  least  squares  or  total^  least  squares  solution 
[20,  p.  222]  of  Eq.  (5.25)  be  and  #2,  then  the  ith 
diagonal  elements  of  and  4>2  represent  the  two  ele¬ 
ments  of  the  ith  parameter  vector,  for  1  ^  i  <  M. 

5.4.  Selection  of  g 

Define  column  vectors,  ef  =  for 

1  <  i  <  M.  Then  X,(g)  =  g^e,-,  for  1  ^  ^  M,  are  the 

eigenvalues  of  (BsBs)~^BsD{g)Bs,  Since  (c^"^s 
for  1  <  Mj  are  distinct  by  assumption,  e-  -  e^  0 

whenever  i  k.  In  order  that  Xi(g),  for  1  <  i  M,  are 
distinct,  it  must  follow  that  g^e^  ¥=  g^e^  or  g^(ej  -  ej^)  ¥=  0, 
for  any  i,  k  with  i  ^  k,  which  implies  that  g*  (the 
superscript  *  denotes  complex  conjugate)  must  not  be 
orthogonal  to  e^  -  e^  for  I  ^  i,  M  and  i  ^  since  Cj 
—  e^j  0  whenever  i  =5^  k.  That  is,  A  in  Eq.  (5.15)  has 
only  distinct  eigenvalues  if  g  does  not  belong  to  the 
M(M  -  l)/4  hyperplanes  expressed  by  g^(e,  -  e^)  =  0 
for  1  ^  i,  k  ^  M  and  i  ¥=  k.  Thus,  if  g  is  selected  ran¬ 
domly,  the  probability  that  A  in  Eq.  (5.15)  has  non- 
distinct  eigenvalues  is,  theoretically,  zero.  (In  compu¬ 
tation,  we  must  consider  some  numerical  margin.)  It 
is  good  to  select  g  randomly.  Then  g'  is  selected  to  be 
orthogonal  to  g,  i.e.,  (g')^g  =  0.  If  e,  -  e^^  is  orthogonal 
to  g*y  then  it  must  be  not  orthogonal  to  (g')*.  Thus, 
Xj(g')  and  X^(g')  become  different  although  Xj(g)  and 
X*(g)  were  the  same. 

Considering  the  fact  that,  if  g  is  selected  randomly, 
the  probability  that  A  in  Eq.  (5.15)  has  non-distinct 
eigenvalues  is  zero  theoretically,  we  can  establish  an¬ 
other  algorithm.  First,  for  L(^2)  randomly  generated 
[^1  ^2  solve  the  Eq.  (5.17)  and  obtain  A's  by  A  = 
BsTig).  Second,  choose  A  which  either  minimizes 

[7,8] 


M 

(5.26) 


or  maximizes 


M 

S5  =  2  a"(e^‘“'S  e^‘^)B5Bfa(c^‘“^S  (5.27) 

i-l 


5.5.  The  Algorithms 


where 

Sxa  =  [SuS2.i*  •  •  •sUa-lSa.Pj-l*  *  •sK-2A-l]'' 

Sib  ~  [82.183,1*  *  *8pj-i.i*  *  •  8^J>j-i83^j-1  *  *  *Sp,-ij>j-i] 

Sla  -  [s'i.iS2,r  *  'Sp,-!,!*  *  *  *  *  •  Spj_i ^.^-2] 


According  to  the  discussion  for  g  in  Section  5.4,  we 
can  set  up  three  algorithms  for  estimating  2-D  param¬ 
eter  vectors.  Algorithm  1  tests  the  magnitude  of  the 
difference  (referred  to,  henceforth,  as  distance)  be¬ 
tween  two  complex-valued  eigenvalues  to  find  if  there 
are  non-distinct  eigenvalues  in  the  matrix  A  of  Eq. 


(5.15).  Algorithm  2  randomizes  the  set  of  eigenvalues 
and  finds  distinct  eigenvalues  by  calculating  the  dis¬ 
tances  between  the  eigenvalues.  Algorithm  3  ran¬ 
domizes  the  set  of  eigenvalues  and  finds  the  set  of  2-D 
parameter  vectors  by  examining  their  fitness  ob¬ 
tained  by  evaluating  the  objective  function  of  Eq. 
(5.26)  or  (5.27). 

In  the  following,  it  is  assumed  that  a  sample  data 
matrix  X  has  been  given  as  in  Eq.  (2.2),  the  number  M 
of  2-D  wavenumbers  is  known,  and  P^{>M)  and 
P2{>M)  have  been  chosen. 

5.5.i.  The  1-D  Based  2-D  PESS  Method  with 
Distance  Test  of  Eigenvalues 

Step  1.  Construct  from  the  sample  data  matrix 
Xby  Eqs.  (3.13)  and  (3.15).  Compute  the  correlation 
matrix  R  by  Eq.  (3.16).  Compute  the  eigendecomposi- 
tion  of  R  to  get  R  =  Partition  B  as  B  =  [Bj^Bg], 

where  is  the  (P1P2)  X  (P1P2  ^  M)  matrix  of  the 
eigenvectors  associated  with  the  smallest  (P1P2  ^  M) 
eigenvalues.  The  sequence  of  eigenvalues  in  A  may  be 
tested  to  estimate  the  number  of  2-D  wavenumbers. 

Step  2.a.  Construct  X^i  and  from  the  sample 
data  matrix  X  by  Eqs,  (4.1),  (4.7),  (4.9),  and  (4.12). 
Compute  the  correlation  matrices  Ri  and  R2  by  Eqs. 
(4.8)  and  (4.13).  Compute  the  eigendecompositions  of 
Ri  to  get  Ri  =  BikiB^,  Partition  Bi  as  Bi  = 
where  is  the  Pj  X  (P^  —  M)  matrix  of  the  eigenvec¬ 
tors  associated  with  the  smallest  (Pj  —  M)  eigen¬ 
values.  Two  univariate  polynomial  equations  in  z,-,  for 
i  —  1,  2,  are  selected  as 

{(B^  + = 

where  1  is  the  column  vector  of  appropriate  order, 
whose  elements  are  all  1. 

Step  2.b.  Construct  two  companion  matrices  D2 
and  Dj  by  Eqs.  (5.8.a),  (5.8.b),  (5.8.c).  Compute 

where  I  is  the  identity  matrix  of  order  Pj, 

G2  =  {B^BsT^B^iD,  ©  /), 

where  I  is  the  identity  matrix  of  order  Pi  and 

G2^{B^Bsr^B^{D2®D,). 

Step  3.  Select  and  g^^^(?^0)  such  that  they 

are  orthogonal.  Also  select  a  threshold  value  5. 

Step  4.  Compute  F  -  +  ^2^^G2  +  with 

g^^^  Compute  the  eigendecomposition  of  F  to  get  F  = 


T*  Diag[Xi  -  •  •  Xa/]  *  P'S  where  the  eigenvector  matrix 
is  [ti-  •  -tA/]. 

Step  5.  Compute  dist,j  =  |  X^  —  X^  | ,  for  1  <  i,  j  <  M, 
and  i  >  j.  Construct  the  clusters  6i  of  eigenvalues  such 
that  \a  —  0\  >  dfoT  aE(Pf^  and  jS E  (Pj  whenever  k¥=j 
and  I  a  —  ^  I  ^  d  for  a,  /S  E  v  Suppose  that  the  num¬ 
ber  of  clusters  is  n. 

Step  6.  If  the  ith  cluster  has  one  member,  keep  tj 
unchanged,  where  is  the  eigenvector  associated  with 
the  eigenvalue  in  the  ith  cluster. 

Step  7.  If  the  ith  cluster  has  m  members,  compute 
Fi  =  (TfT,)-^Tf(^i2>Gi  +  G2  +  G3)T,  with 

where  Tj  is  the  matrix  of  the  eigenvectors  associated 
with  the  eigenvalues  in  the  ith  cluster.  Compute  the 
eigenvector  matrix  Ui  of  Pj  and  replace  Tj  by  Tj  L/’j. 

Step  8.  Compute  the  estimate  A  of  A  as  A  =  PsT. 
Construct  Si^,  Sn,,  Sza,  ^26  by  Eq.  (5.25).  Solve  the 
linear  systems  via  the  method  of  least  squares  to  get 

-Ji  =  (5f„5 

and 

^2  =  (SgS2X^SgS2,. 

(^2)u)}i^x  is  the  set  of  2-D  signal  zeros  which 
give  the  estimate  of  2-D  wavenumbers. 

5.5.2.  The  1-D  Based  2-D  PESS  Method  with 

Distance  Test  of  Randomized  Eigenvalues 

Step  1,  2.  Same  as  that  of  the  algorithm  in  Sub¬ 
section  5.5.1. 

Step  3.  Initialize  outer  loop;  m  =  1,  G^f^  =  Gi,  G^^^  = 
G2,  G^^^  =  G3  and  D  =  I,  where  I  is  the  identity  matrix 
of  order  M.  Choose  Tq*  and  rjj,  for  1  <  f  <  3,  in  ^j  = 
with  Toj  ^  Tj  <  Tij. 

Step  4.  Initialize  inner  loop;  n  =  0,  d _ =  -1. 

Step  5.  Generate  a  complex  random  vector  (gi,  g2y 
^3).  Compute  P  =  +  ^2G2'”^  +  g^Gz'^-  Compute 

the  eigendecomposition  F  -  Diag[Xi  •  •  •  X^^]  •  IPh 

Step  6.  Compute 

d  =  max  min  |  Xj  ~  X*  | 

l<i<M 

and  denote  the  maximizing  argument  of  i  by  If  d  > 
dniax  then  set  d^  =  d,  =  q  and  T  =  L7. 

Step  7.  Compute  n  =  n  +  l.  Ifn<L  then  go  to 
step  5. 


Step  8.  Rearrange  the  columns  of  T by  T  =  [tiT2], 
where  ti  is  the  column  of  T  and  Tj  is  the  matrix 

of  all  the  remaining  columns.  Replace  D  by 
D-Diag[/„.i,  T]. 

Step  9,  Compute  m  =  Ifm<M+l  then 

compute  =  (TfT2)"^T^G,T2  for  1  ^  i  <  3,  and  go 
to  Step  4. 

Step  10.  Same  as  Step  8  of  the  algorithm  in  Sub¬ 
section  5.5.1. 

5.5,3.  The  1-D  Based  2~D  PESS  Method  with 
Fitness  Test  of  2-D  Parameter  Vectors 

Step  1,  2.  Same  as  that  of  the  algorithm  in  Sub¬ 
section  5.5.1, 

Step  3.  Initialize  the  loop;  m  -  0,  d^^^  =  — 1. 
Choose  and  for  1  ^  i  <  3. 

Step  4.  Generate  a  complex  random  vector  (gi, 
g^)  where  gi  =  and  r,o  ^  r,-  <  r^^.  Compute  F  =  g^G^ 
82^2  ■*"  8z^3*  Compute  the  eigendecomposition  of  F 
byF=T-Diag[X,-..X^].T-\ 

Step  5.  Compute  A  =  BsT. 

Step  6.  Construct  52^  and  ^25  by  Eq. 

(5.25).  Solve  the  linear  system  via  the  least  squares 
method  by 

and 

Denote  the  set  of  2-D  signal  zeros  (^2)u)}i!ii 

by  c^. 

Step  7.  Compute  d  =  Sill  where  a,-  ~ 

a(('^i)»i  (^2)«)>  where  a(  • )  is  defined  by  Eq.  (3.5).  If  d 
>  then  set  G  =  and  =  d. 

Step  8.  Compute  m  =  m  +  l.  Ifm<L  +  l  then  go 
to  step  4. 

Step  9.  Each  2-D  signal  zero  in  the  set  &  gives  the 
estimate  of  a  2-D  wavenumber. 

^HE  PAIRING-PESS  METHOD 

The  algorithm  presented  in  the  previous  section 
can  be  used  as  an  algorithm  for  pairing  separate  1-D 
parameters  estimated  by  1-D  based  2-D  methods.  As¬ 
sume  that  we  have  estimated  1-D  signal  zeros, 
and  and  we- know  the  number  M 


of  2-D  sinusoids.  Using  the  signal  zeros  estimated, 
construct  the  polynomial  equations  by 

2  ^  (zi  -  ri  (zi  -  (6.1.a) 

At-0  A-1 

Pi 

2  (Z2  -  1)“^  n  (Z2  -  (6.1.b) 

A-0  *-l 

where  dj  =  0  if  M  =  M,,  or  dj  =  1  if  M  >  Mj,  for  i~  1,2. 
When  M  >  M,,  some  of  e*'"'**,  for  1  ^  are 

multiple  order  zeros.  The  case  of  multiple  1-D  signal 
zeros  is  accounted  for  by  the  2-D  Vandermonde  sys¬ 
tem.  Suppose  that  data  matrix  X  in  Eq.  (2.2)  for  the 
case  when  M  =  2  in  Eq.  (2.1)  was  available.  Suppose 
that  the  estimated  1-D  signal  zero  sets  are  {e^'^}  and 
Here  Mi  =  1,  and  Mg  =  1  and  using  these  two 
sets  of  1-D  signal  zeros,  Eq.  (6.1.a)  and  (6.1.b)  are 
written  as 

zf‘“^-(zi  -  l)(zi  -  e-''*)  =  0,  (6.2.a) 

zf^-^-(z2-l)(z2-e^-")  =0.  (6.2.b) 

When  we  construct  the  2-D  Vandermonde  system, 
using  the  two  1-D  polynomial  equations  in  Eq.  (6.2.a) 
and  (6.2.b),  described  by  Eq,  (5.13),  a(l,  1),  a(l,  e^’^), 
a(e^'\  1),  and  a(e-''\  e^'^)  are  the  eigenvectors  of  D(g). 
Theorem  5.1  enables  us  to  extract  the  subset  of  eigen¬ 
vectors  associated  with  the  two  2-D  wavenumbers. 
The  remaining  eigenvectors  account  for  other  incor¬ 
rectly  paired  2-D  zeros  and  spurious  2-D  zeros.  When 
the  algorithm  given  in  Section  5  is  used  as  a  general 
pairing  algorithm,  Eq.  (6.1.a)  and  (6.1.b)  are  substi¬ 
tuted  for  Eq.  (5,4)  and  (5.7).  All  the  remaining  proce¬ 
dure  is  used  with  no  change. 

6.1.  The  Algorithms 

In  the  following,  it  is  assumed  that  a  sample  data 
matrix  X has  been  given  as  in  Eq.  (2.2),  the  number  M 
of  2-D  wavenumbers  is  known,  and  Pi{>M)  and 
P2{>M)  have  been  chosen. 

6.LL  The  Pairing-PESS  Method  with  Distance 
Test  of  Eigenvalues 

The  steps  are  the  same  as  those  of  the  algorithm  in 
Subsection  5.5,1,  except  that  Step  2. a  is  changed  to 
the  following. 

Step  2.a.  Suppose  that  the  estimates  of  the  sepa¬ 
rate  frequencies  are  given  by  the  sets  and 

{ojgjjili,  where  Mj  <  M  and  Mg  <  M.  Two  univariate 
polynomial  equations  are  given  by  Eqs.  (6.1. a)  and 
(6,l.b), 


6.1.2.  The  Pairing'PESS  Method  with  Distance 
Test  of  Randomized  Eigenvalues 

The  steps  are  the  same  as  those  of  the  algorithm  in 
Subsection  5.5.2,  except  that  Step  2.a  is  changed  to 
the  Step  2.a  of  the  algorithm  in  Subsection  6.1.1. 

6.1.3,  The  Pairing-PESS  Method  with  Fitness 
Test  of  2-D  Parameter  Vectors 

The  steps  are  the  same  as  those  of  the  algorithm  in 
Subsection  5.5.3,  except  that  Step  2.a  is  changed  to 
the  Step  2.a  of  the  algorithm  in  Subsection  6.1.1. 

7.  DERIVATION  OF  THE  DIRECT  2-D  PESS 
METHOD 


In  this  section,  it  is  shown  that  the  2-D  wavenum¬ 
bers,  (wip  ai2i),  for  1  ^  i  <  M,  are  directly  obtainable  by 
using  the  structure  of  the  vector  a(e‘^"S  and  the 
estimate  of  the  range  space  of  the  matrix  A  defined  in 
Eq.  (3.6). 

7.1.  A  Basic  Matrix  Equation 

For  convenience,  the  following  convention  for  in¬ 
dexing  the  elements  of  a  {P^  •  Pi)  X  1  column  vector  is 
introduced.  The  {Pi  •  {j  —  1)  +  i}th  element  of  the 
vector,  where  1  ^  7  <  P2  and  1  <  f  <  Pi,  will  be  called 
the  (7,  i)th  element.  Thus,  the  (7,  0th  element  of  the 
(P2'Pi)  X  1  column  vector  a(zi,  Zg)*  after  replacing 
by  Zi  and  by  Zj,  is 

=  z{“^-zr^  (7.1) 

The  elements  of  a(zi,  Zg)  are  then  shown  explicitly 
below. 

a(zi,  22) 

=  [aiaOi^i-  •  -ai^j-  •  (7.2) 

From  Eq.  (7.1),  it  follows  that 

(§2^2  El^l)  “  ^  §l^j.i+l9  C7*3) 

where  1  <  i  <  Pi  1, 1  <7  <  P2  ~  1,  and^i  and^a  are 
arbitrary  complex  numbers.  Let  e^  ^  be  the  (7,  0th  col¬ 
umn  of  the  identity  matrix  of  order  P2  •  Pi,  i.e.,  the 
(P2  •  Pi)  X  1  column  vector  whose  (7,  0th  element  is 
one  and  all  the  remaining  elements  are  zero.  Then, 
Eq.  (7.1)  can  be  expressed  as 

(7.4) 


Eq.  (7.3)  can  then  be  reformulated  as 

ig^Zi  +  giZi)  •  e^a(zi,  Z2) 

“  ^2)  22).  (7.5) 

where  1  <  i  <  Pi  —  1  and  1  ^7  <  P2  —  1.  Now,  let  e^^j,  for 
1  ^  i  <  Pi  1  and  1  7  <  P2  1,  be  stacked  as  the 

matrix, 

c  =  [euCu-  •  -eij.,-!-  •  • 

(7.6) 

which  has  (P1P2)  rows  and  (Pi  —  1)(P2  ~  1)  columns. 
Stacking  both  sides  of  Eq.  (7.5),  it  follows  that 

(^2^2  +  Si^i)  •  C^a(zi,  Z2)  =  D^igu  ^2)a(2i,  Z2),  (7.7) 

where 

Digiy  gi)  =  g2D2  +  giDu  (7.8) 

D2  ^  [®2,1®2.2*  *  *®2.Pi-I*  *  *  '  *  *®P2hPi-iT  (7.9) 

and 

Di  =  [eu®u*  •  •  •  ap2-u®P2-i.3 ’  ‘  •®P2-i.pJ* 

(7.10) 

The  matrix  equation  (7.7)  will  be  used  to  develop  the 
new  algorithm.  Eq.  (7.7)  holds  for  an  arbitrary  2-tuple 
(zi,  22),  as  is  obvious  from  the  steps  in  its  derivation. 
Substituting  (e^"‘*,  e*'"^)  for  (zi,  22)  in  Eq.  (7.7)  yields 

(7.11) 

where  1  ^  i  ^  M  and 

\igi9  g2)  =  g2e^^  +  (7.12) 

Through  the  definition  of  the  matrix  A  extracted 
from  Eq.  (3.6),  Eq.  (7.11)  leads  to 

D^(Mx9  g2)A  =  C^AA(^i,  ^2),  (7.13) 

where  A(^i,  is  a  diagonal  matrix  whose  ith  diagonal 
element  is  X,(^i,  g^f)  and  whose  order  is  M,  Eq.  (7.13) 
can  be  solved  for  A  using  the  property  of  signal  (zeros) 
selectivity  in  the  signal  subspace,  as  described  in  the 
next  section. 


°j.i  =  22)- 


7.2.  Estimation  of  the  Matrix  A 

Substituting  Eq.  (3.19,b)  for  A  in  Eq.  (7.13),  the 
problem  of  solving  for  A  is  reduced  to  the  problem  of 
obtaining  the  M  X  M  nonsingular  matrix  T  in 

g2)Bs}  •  T  =  {C^Bs)  •  T.  g^)  (7.14) 

Using  Eqs.  (3.6)  and  (7.6),  it  can  be  shown  that 

C^A  =  (7.15) 

where 

aXe^"“,  =  a^(e-'‘^0  ®  ai(e^“'0, 

and 

ai(e^"'0  =  - 

From  the  discussion  in  Section  3.3  and  [4],  it  is  known 
that  C^A  is  of  full  rank  if  and  only  if 

(P,  -  1)  •  (P2  -  1)  ^  M.  (7.16) 

Condition  (7.16)  is  always  satisfied  since  it  was  as¬ 
sumed  that  Pi  >  M  and  P2  >  M.  Since  A  and  P5  span 
the  same  signal  subspace,  C^Bs  is  also  of  full  rank. 
Thus,  Eq.  (7.14)  implies  that 

F{g,,g^)  =  T^A(g,,g,)-T-\  (7.17) 

where 

F{g^,  g2)  =  (B^CC^Bsr^B^CD^igu  ^2)^5.  (7.18) 

The  M  X  M  non-singular  matrix  T  and  the  diagonal 
matrix  A(^i,  ^2)  be  obtained  by  eigendecomposing 
the  matrix  Figi,  ^2).  Eq.  (7.17)  implies  that  the  matrix 
F  is  diagonalizable  (non-defective)  [20,  p.  338]. 

Lemma  7.1.  With  a  distinct  eigenvalue  is  asso¬ 
ciated  one  eigenvector  which  is  unique  up  to  scalar 
multiplication. 

Proof.  F  E  can  be  Schur-decomposed  as 

Q^PQ  ==  T,  where  Q  E  is  a  unitary  matrix  and  T 
E  is  a  upper  triangular  matrix  [20,  pp.  335]. 
Since  det(P  -  XI)  =  det(QTQ^  -  X/)  =  det(T  -  XT), 
the  diagonal  elements  X,-,  for  l^i^  M,  are  the  eigen¬ 
values  of  F.  Consider  an  eigenvalue  X^.  Its  associated 
eigenvector  x,*  satisfies  (P  —  X^J^x^  =  0  or  (T  “  X(/)Q^x, 
=  0.  If  Xj  is  distinct,  T  —  \I  is  of  rank  M  —  1.  Thus, 
Q^Xj  is  unique  up  to  scalar  multiplication,  and  so  is 
x,.H 


When  the  eigenvalues  of  the  matrix  P(^i,  ^2)  ^^e 
distinct  (P(^i,  ^2)  is  nonsingular  if  and  only  if  any 
eigenvalue  is  not  zero.),  the  columns  of  the  eigenvec¬ 
tor  matrix  T  are  unique  up  to  scalar  multiplication, 
and  so  the  estimate  of  A  is  given  by 

A  =  P5TA,  (7,19) 

where  A  is  a  diagonal  matrix  whose  elements  are  de¬ 
termined  to  scale  each  column  of  P5T  such  that  the 
first  element  is  1.  The  estimate  of  2-D  wavenumbers 
are  computed  from  Eq.  (7.19)  using  Eqs.  (5.23),  (5.24), 
and  (5.25). 

7.3.  The  Algorithms 

In  the  following,  it  is  assumed  that  a  sample  data 
matrix  X  has  been  given  as  in  Eq.  (2.2),  the  number  M 
of  2-D  wavenumbers  is  known,  and  Pi{>M)  and 
P2(>M)  have  been  chosen. 

7.3.1.  The  Direct  2-D  PESS  Method  with 
Distance  Test  of  Eigenvalues 

The  steps  are  the  same  as  those  of  the  algorithm  in 
Subsection  5.5.1,  except  that  ^3  and  G3  vanish  in  all 
the  steps  and  the  Steps  2.a  and  2.b  are  changed  to  the 
following. 

Step  2.  D2  and  are  constructed  by  Eqs.  (7.9) 
and  (7.10),  and  the  following  2  matrices  are  com¬ 
puted. 

Gi  =  {B^CC^BsT^B^QD^Bs 

and 

G2  =  (B^CC^BsT^B^CD^Bs. 

7.3.2.  The  Direct  2-D  PESS  Method  with 
Distance  Test  of  Randomized  Eigenvalues 

The  steps  are  the  same  as  those  of  the  algorithm  in 
Subsection  5.5.2,  except  that  ^3  and  G3  vanish  in  all 
the  steps  and  the  Steps  2.a  and  2.b  are  changed  to  the 
step  2  of  the  algorithm  in  Subsection  7,3.1. 

7.3.3.  The  Direct  2-D  PESS  Method  with 
Fitness  Test  of  2-D  Parameter  Vectors 

The  steps  are  the  same  as  those  of  the  algorithm  in 
Subsection  5.5.3,  except  that  ^3  and  G3  vanish  in  all 
the  steps  and  the  Steps  2.a  and  2.b  are  changed  to  the 
Step  2  of  the  algorithm  in  Subsection  7.3.1. 


8.  SIMULATION  RESULTS 


as 


Step  2.b.  The  initial  G„  and  G3  were  computed 


In  this  section,  examples  are  presented  to  illustrate 
the  PESS  methods  developed  above,  and  Monte  Carlo 
simulation  results  are  provided  to  show  the  perfor¬ 
mance  of  the  methods.  In  the  simulation,  the  PESS 
method  and  the  MEMP  method  [8]  were  applied  to 
the  same  data  set  for  accuracy  comparison. 

Each  20  X  20  sample  data  matrix  was  generated  by 
the  equation 

3 

x(m,  n)  =  2  +  v(m,  n), 

A-X 

for  0  ^  m,  n  ^  19,  (8.1) 

where  v{m,  n)  is  the  complex  white  noise. 

8.1.  Examples 

The  3  wavenumber  parameters  were  selected  as 

(/ii,/2i)  =  (0.26,0.24)  (8.2) 

(/i2,/22)  =  (0.24,0.24)  (8.3) 

(/i3,/23)  =  (0.24,  0.26).  (8.4) 

The  3  algorithms  processed  a  sample  data  matrix  X 
which  was  generated  with  SNR  of  20  dB.  SNR  is  de¬ 
fined  by 

SNR  =  10  logio(l/<^),  (8.5) 

where  is  the  noise  power.  The  4X4  array  was  se¬ 
lected,  which  means  that  —  P2  ~  The  number  L 
of  iteration  was  set  to  2  and  Vqi  =  1  and  =  10,  for  1  < 
i  <  3,  were  selected  in  Sections  8.1.2  and  8.1.3. 

8.1.1.  The  1-D  Based  2-D  PESS  Method  with 
Distance  Test  of  Eigenvalues 

Step  1.  The  list  of  data  matrix  X,  R,  Pi,  P2»  By  ^y 
Pi,  Ai,  P2,  and  Ag  were  obtained  and  is  available  from 
the  authors  as  an  Appendix  to  this  manuscript. 

Step  2.a.  Two  univariate  polynomial  equations  in 
Zi,  for  i  =  1,  2,  were  computed  as 

0.0  +  (-0.0204075  +  ;0.1130556)Z|  +  (-0.2097420  +  j0.0769447)rJ 
+(-0,2097420  -  i0.0769447)*f  +  (-0.0204075  -  j0.U30556)zJ  =  0.0 

and 


0.0236260  +  ;0.9975650  0.0321896  -  jO.0225824  -0.0000030  -  ;0.0015337\ 

0.0686103  +  i0.0446103  -0.0180299  +  >1.0018485  -0.0025797  +  j0.0040438 
0.0155744  +  >0.6278006  -0.3369688  -  >0.5409832  0.0179116  +  >0.9976731  ) 

02  = 

0.0177350 +  >0.9998431  -0.0400412 +  >0.0150251  0.0001062  -  >0.001 1839  \ 

-0.0760086  -  >0.0483654  -0.0297496  +  >0.9890730  0.0025833  -  >0.0035922 
0.0096372  +  >0.6282084  0.3278028  +  >0.5331685  0.0181662  +  >0.9980672/ 


and 


(73  = 

- 1.0044686  +  >0.0348798  0.0043723  -  >0.0093423  0.0027056  -  >0.0000879  \ 

0.0038265  -  >0.0068570  -0.9936933  -  >0.0474929  -0.0003874  -  >0.0000307  . 

-1.2538954  +  jO.  1007968  0.0082100  -  >0.0070122  -0.9974769  +  >0.0359509  / 


Step  3.  g  =  [0  0 1],  g'  =  [1  -0.5  0]  and  5  =  0.07  were 
selected. 

Step  4.  F  =  giGi  -I-  giGz  +  gs,G^  was  computed  as 


F  = 

(-1.0044686  +  >0.0348798  0.0043723  -  >0.0093423  0.0027056  -  >0.0000879  \ 

0.0038265  -  >0.0068570  -0.9936933  -  >0.0474929  -0.0003874  -  >0.0000307 

-1.2538954  +  >0.1007968  0.0082100  -  >0.0070122  -0.9974769  +  >0.0359509  / 

The  eigenvalue  and  the  eigenvector  matrices  of  F 
were  obtained  as 

A  =  ato9(- 1.0048345  -  >0.0229922  -  0.9980068  +  >0.0940336  -  0.9927975  -  >0.0477037] 


and 


(0.0021124 +  >0.0471462  0.0041338 ->0.0459490  -0.0075787  +  >0.0682349  \ 

0.0054370 ->0.0001844  -0.0015099  +  >0.0048979  -0.0158513  +  >0.3115187  j  . 
l.OOOOOOO  +  >0.0000000  1.0000000  +  >0.0000000  1.0000000  +  >0.0000000  / 

Step  5.  The  distances  between  the  eigenvalues 
were  computed  as  in  the  distance  matrix 

(0.0000000  0.1172248  0.0274872  \ 

0.1172248  0.0000000  0.1418330  . 

0.0274872  0.1418330  O.OOOOOOO/ 

The  clustering  procedure  gave  cluster  1:  {2}  and  clus¬ 
ter  2:  {1,3}. 

Step  6.  Consider  cluster  1.  Keep  the  second  col¬ 
umn  of  T. 

Step  7.  Consider  cluster  2.  The  following  F'  was 
computed  with  g'  =  [1  —0.5  0]. 

p,  _  f  0.0430450  +  >0.5375658  0.0366425  +  >0.0193354 

^  V  0.2040451 ->0.1839002  -0.0454788  +  >0.4761306/ ‘ 


The  eigendecomposition  F'  =  T'N(T')  ^ 
puted  as 


0.0  +  (0.3638633  -  >0.4630003)^,  +  (0.6093172  +  >0.2134817)z| 
+(0.6093172  -  >0.2134817)x^  +  (0.3638633  +  >0.4630003)i^  =  0.0. 


was  com- 


A'  =  Dta  j[-0.1 109810 +  ;0.50718S40. 10854 73 +  J0.5065  no] 


and 


-0.2528226  -  >0.0756661  0.3424790  +  >0.4575574 

1.0000000  +  j 0.0000000  1.0000000  +  >0.0000000  ;  • 


Step  8.  The  desired  eigenvector  matrix  was  ob¬ 
tained  as 


/ -0.0317045  +  j 0.06704 75  0.0041338  ^  >0.0459490  0.0091916  +  >0.0304683\ 

-0.1425295 +  >0.0992513  -0.0015099  +  >0.0048979  0.0330160  -  >0.0777440 
\  1.3424790  +  >0.4575574  1.0000000  +  >0.0000000  0.7471774  -  >0.0756661  J 


The  estimate  of  A  was  computed  as 


A  = 

0  3513742  +  >0.0823592  0.2483042  +  >0.0218955  0.1825521  -  >0.0333121  \ 

-0  0606784  +  >0.3559272  -0.0092024  +  >0.2492167  0.0184874  +  >0.1869526 

-0  3587124  -  >0.0390216  -0.2493843  +  >0.0033865  -0.1901246  +  >0.0043823 

0  0153205 ->0.3609187  -0.0156432  -  >0.2488627  0.0108940  -  >0.1925282 

-0.1142487  +  >0.3388163  -0.0107822  +  >0.2498093  0.0452373  +  >0.1799749 

-0.3464285  -  >0.0922990  -0.2499995  +  >0.0016945  -0.1853295  +  >0.0305067 
0  0718142  -  >0.3515703  -0.0142784  -  >0.2497213  -0.0168328  -  >0.1891642 

0.3557070 +  >0.0490836  0.2486869  -  >0.0266547  0.1924566  -  >0.0018425 

-0.3256860  -  >0.1427773  -0.2503781  +  >0.0002922  -0.1771500  +  >0.0557414 
0.1216560 ->0.3349905  -0.0127796  -  >0.2502443  -0.0414145  -  >0.1831561 

0.3420828  +  j0.1018271  0.2493742  -  >0.0253760  0.1880830  -  >0.0280891 

-0  0794458  +  >0.3480729  0.0377740  +  >0.2480116  0.0131399  +  >0.1921169 

0.1726683 ->0.3085556  -0.0111399 ->0.2505019  -0.0674049  -  >0.1734022 

0.3199914 +  >0.1519383  0.2498546  -  >0.0235370  0.1802734  -  >0.0533663 

-0.1325859  +  jO.3296704  0.0361808  +  j0.2485306  0.0400952  +  >0.1859581 

,  -0.3379245  -  >0.1122691  -0.2464574  +  >0.0488280  -0.1908892  +  >0.0256973/ 


From  A,  the  following  two  matrices  were  obtained. 


♦i  = 

/  0.0662453  +  >0.9975544  0.0008960  -  >0.0005833  -0.0012975  -  >0.0003834  \ 

1  0.0022932  +  >0.0024737  0.0493609  +  >0.9988683  -0.0021381  +  >0.0109805  i 

\  -0.0099251  -  >0.0038398  -0.0010368  +  >0.0006666  -0.0728221  +  >0.9969979/ 

/  -0.0835589  +  >0.9959845  -0.0023614  +  >0.0024702  0.0025242  -  j0.0003328  \ 

f  -0.0182159  -  >0.0148995  0.0472167  +  >0.9988186  -0.0016406  -  >0.0004648  j 

\  0  0119660  +  >0.0059019  0.0012700  -  >0.0000055  0.0606975  +  >0.9979546  / 


The  desired  2-D  zeros  were  extracted  from  the  above 
two  matrices  as 


(0.9997516  •  «p(;27r0.2394464),  0.9994834  ■  «p(j2t0.2633212)), 
(1.0000872  •«2p(>2jr0.2421415),  0.9999340- exp(>2T0.2424819)), 
(0.9996539  •  exp(>2jr0.2616043),  0.9997987  •  exp(>2T0.2403318)). 


8.1.2.  The  Pairing-PESS  Method  with  Distance 
Test  of  Randomized  Eigenvalues 
It  was  assumed  that  two  sets  of  1-D  wavenumber 
parameters  were  estimated  as 

{0.2595620,  0.2396049}, 

{0.2567192,  0.2431879}. 


Two  univariate  polynomial  equations  in  Zi,  for  i  —  1, 2, 
were  obtained  as 

0.0  +  (0.9999863  -  >0.0052345)zi  +  (-0.9947621  +  >2.0012980)z? 


+(-1.0052242  -  >l.9960635)x?  +  l.OzJ  =s  0 


and 

0.0  +  (0.9999998  -  >0.0005840)rj  +  (-0.9994164  +  >l.9987771)z3’ 

+(-1.0005834  -  >1.9981931)z|  +  l.Oz^  =  0. 

Iteration  1.  The  vector  g  was  randomly  generated 
as 

g  =  {8.2824010  +  >2.7776125  1.7048545  -  >4.0095459  0.3119416  -  >0.9788936]'. 

The  maxmin  distance  was  computed  as  1.0368955. 
Another  g  was  randomly  generated  as 

g  =  [-0.1596835  +  >2.1919652  6.8230877  +  >7.0010181  -  6.1613831  -  >8.6783894]'. 

The  maxmin  distance  was  computed  as  1.0061741. 
The  first  trial  was  determined  to  be  better.  The  ma¬ 
trix  D  was  obtained  as 

{  0  0068896  +  >0.0257567  0.0072275  +  >0.0205755  0.0029841  -  >0.0483550  \ 

0  0476577  -  >0.0772910  -0.0524261  +  >0.0925424  -0.0107985  +  >0.0038552  . 

\  1.0000000  +  >0.0000000  1.0000000  +  >0.0000000  1.0000000  +  >0.0000000  j 

Iteration  2.  The  vector  g  was  randomly  generated 
as 

g  -  [5.4793580  +  >6.2770135  -  1.5739585  +  >0.5342874  2.0556196  -  >3.2626202]'. 

The  maxmin  distance  was  computed  as  1.6158669. 
Another  g  was  randomly  generated  as 

g  ==  [5.8321917+>7.7009918  -0.2320668+>3. 1577547  ^3  =  -4.355379l->9.0153360]'. 

The  maxmin  distance  was  computed  as  1.8360912. 
The  matrix  D  was  updated  as 

/  0.0068896  +  >0.0257567  0-0028700  -  >0.0482168  -0.0066433  -  >0.0632031  \ 

(  0.0476577 ->0.0772910  -0.0116724  +  >0.0038789  0.0128859  -  >0.0791381  . 

\  l.OOOOOOO  +  >0.0000000  1.0042443  +  >0.0070381  0.2113152  +  >0.1908670  j 

The  desired  2-D  zeros  were  extracted  from  the  above 
two  matrices  as 

(1.0003608-  exp(>2T0.2615938),  1.0001100  •  exp(> 2x0.2402652)), 

(1.0002271  ■  €ip{>2x0.2419204),  0.9929687-  cxp<>2x0.2431732)), 

(0.9988808  •  exp(>2x0.239662I),  1.0055987  •  exp(>2x0.2625731)). 

8.1.3.  The  Direct  2-D  PESS  Method  with 
Fitness  Test  of  2-D  Parameter  Vectors 

Trial  1.  The  following  were  generated  randomly. 

si  =  8.2824010  +  >2.7776125  aad  s2  =  1.7048545  -  >4.0095459. 

The  eigenvector  matrix  was  computed  as 


-0.0022579  +  ;0.0378996  -0.0061035  +  ;0.032090  -0.0024479  -  j0.0504617 
0.0563264  -  j0.0876339  -0.0741348  +  j0.087846  0.0008896  -  >0.0042041 
1.0000000 +  >0.0000000  1.0000000  + >0.000000  1.0000000  +  >0.000)000 


With  this  estimate  of  T,  the  signal  zeros  were  com¬ 
puted  as 

(-0.0729085  +  >0.9971211,  0.0611014  +  >0.9980343) 

(0.0666534  +  >0.9975640,  -0.0843249  +  >0.9953308) 

(0.0490392  +  >0.9987356,  0.0475788  +  >0.9993925), 

and  the  criterion  value  was  evaluated  as  47.9998695. 
Trial  2.  The  following  were  generated  randomly. 

j;(l)  =  -0.1596835  +  >2.1919652  and  5(2)  =  6.8230877  +  >7.0010181. 

The  eigenvector  matrix  was  computed  as 

r  = 

/  -0.0075360  +  >0.0414075  -0.0001785  -  >0.033024  -0.0139303  +  >0.0298457\ 

{  -0.0666335 +  >0.1067764  0.0090331  -  >0.015174  0.0623083  -  >0.0678467  1. 

\  1.0000000  +  >0.0000000  1.0000000 +  >0.000000  1.0000000  +  >0.0000000  / 

With  this  estimate  of  T,  the  signal  zeros  were  com¬ 
puted  as 

(0.0664314  +  >0.9976989,  -0.0845125  +  >0.9959812) 

(0.0513499  +  >0.9941209,  0.0470659  +  >0.9994206) 

(-0.0749972 +  >1.0016009,  0.0618018  +  >0,9973559). 

and  the  criterion  value  was  evaluated  as  47.9998661. 

The  result  of  the  first  trial  was  determined  to  be 
best.  The  desired  2-D  zeros  were  extracted  from  the 
above  two  matrices  as 

(0.9997830-  exp(>2jr0.26l6166),  0.9999029  ■  exp(> 2x0.2402684)), 

(0.9997883  •  erp(>2x0.23938l7),  0.9988964  •  exp<>2x0.2634516)), 

(0.9999388  •  ejp(>2x0.2421916),  1.0005244  -  exp(>2x0.2424287)). 


FIG.  8.2.  The  pairing-PESS;  20  dB;  200  runs;  Pi  ^  P2  =  4. 

8.2.  Monte  Carlo  Run 

The  simulation  results  for  the  1-D  based  2-D  PESS 
method,  the  Pairing-PESS  method,  and  the  Direct 
2-D  PESS  method  are  presented.  The  procedure  for 
avoiding  the  problem  of  repeated  eigenvalues  was  se¬ 
lected  to  provide  the  fitness  test  of  2-D  parameter 
vectors.  Comparison  of  the  PESS  methods  with  the 
MEMP  method  was  also  made.  Two  hundred  sample 
data  matrices  were  generated  and  the  three  PESS 
methods  and  the  MEMP  method  were  applied  to  esti¬ 
mate  the  2-D  wavenumbers.  The  comparisons  be¬ 
tween  the  four  methods  were  made  for  the  wavenum¬ 
ber  parameters  in  Eqs.  (8.2),  (8.3),  and  (8.4).  (The 
results  are  shown  in  Figs.  8.1  to  8.4.) 

Subsequently,  the  Direct  2-D  PESS  method  and 
the  MEMP  method  were  compared  for  wavenumber 
parameters  specified  by 

ifu,  /21)  =  (0-26,  0.24)  (8.6) 

(/i2,/22)  =  (0.24,0.24)  (8.7) 

(/i3, /23)  =  (0.25,  0.26).  (8.8) 
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FIG.  8.1.  The  1-D  based  2-D  PESS;  20  dB;  200  runs;  P,  =  Pj  =  5. 
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FIG.  8.3.  The  direct  2-D  PESS;  200  runs;  P,  =  P,  =  4. 


FIG.  8.4.  The  MEMP;  20  dB;  200  runs;  Px^  4. 


FIG.  8.6.  The  direct  2-D  PESS;  20  dB;  200  runs;  Pi  =  P2  =  5- 


For  the  preceding  set  the  parameters  and  P2  were 
each  selected  to  be  5  to  improve  resolution.  The  re¬ 
sults  for  this  case  are  shown  in  Figs.  8.5  and  8.6. 

In  the  simulations,  =  1  and  r^i  =  10,  for  1  <  i  <  3, 
were  selected.  The  number  L  of  iterations  was  chosen 
to  be  2  for  the  PESS  algorithms.  In  the  simulation  of 
the  Pairing-PESS  method,  it  was  assumed  that  the 
estimates  of  wavenumbers,  obtained  by  other  method, 
include  an  estimation  error  with  a  standard  deviation 
of  0.003.  In  the  pairing  procedure  by  the  MEMP 
method,  the  criterion  of  Eq.  (5.27)  was  used  as  recom¬ 
mended  by  Hua  [8]. 

Figures  8.1  to  8.4  show  that  the  PESS  methods 
worked  as  well  as  the  MEMP  method.  The  Direct  2-D 
PESS  method  outperformed  the  MEMP  method  for 
the  2-D  wavenumber  parameters  in  Eqs.  (8.6)  to  (8.8), 
as  shown  in  Figs.  8.5  and  8.6. 

9.  CONCLUSIONS 


The  methods  proposed  here  for  estimation  of  2-D 
wavenumbers  of  a  finite  number  of  sinusoids  in  noise 
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FIG.  8.5.  The  MEMP;  20  dB;  200  runs;  Pj  =  P2  =  5. 


are  computationally  efficient  and  accurate.  The 
knowledge  of  the  number  of  sinusoids  or  its  estimate 
is  needed  but  this  constraint  can  be  relaxed  as  in  the 
1-D  case  [3].  The  procedure,  in  principle,  is  generaiiz- 
able  to  dimensions  higher  than  two.  Though  the  com¬ 
putations  are  likely  to  increase  substantially  with  in¬ 
crease  in  the  number  of  dimensions,  other  1-D  based 
methods  which  use  a  separate  pairing  procedure  as 
well  as  other  direct  methods  are  usually  more  compu¬ 
tation-intensive. 
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ABSTRACT 

A  theoretical  analysis  is  proposed  to  evaluate  the  ro¬ 
bustness  of  the  total  least  square  (TLS)  algorithm  for 
image  reconstruction  from  a  sequence  of  imdersampled 
noisy  and  blurred  frames  in  the  wavenumber  domain. 
Since  the  data  samples  in  the  wavenumber  domain  are 
complex- valued,  the  results  from  the  TLS  theory  de¬ 
veloped  for  the  case  of  real-valued  data  are  adapted  to 
I  the  situation  at  hand.  It  is  shown  that  the  image  qual- 
l  ity  can  be  improved  as  more  frames  become  available. 
I  In  the  case  of  blurred  frames,  higher  resolution  images 
I  may  be  reconstructed  using  the  TLS  algorithm  with 
post-deblurring.  Finally,  computer  simulation  results 
are  provided  to  demonstrate  the  robustness  of  the  TLS 
algorithm  for  image  reconstruction. 

1.  INTRODUCTION 

Recently,  considerable  attention  has  been  devoted  to 
image  sequence  processing.  Many  applications  of  im¬ 
age  sequence  processing  may  be  found  in  (both  military 
and  industrial)  surveillance,  medical  imaging,  and  agri¬ 
culture.  One  of  the  efforts  in  image  sequence  processing 
is  directed  at  the  improving  of  resolution  of  the  recon¬ 
structed  image  from  a  sequence  of  snapshots  taken  over 
a  fixed  object.  For  example,  LANDSAT  follows  a  cir¬ 
cular,  near-polar  orbit  at  a  nominal  altitude  of  700km 
and  it  flies  over  the  earth  surface  between  latitudes  82® 
north  and  82®  south  once  every  16  days.  Thus,  multi¬ 
ple  images  of  the  same  area  are  available.  The  spatial 
resolution  of  an  image  is  often  determined  by  imaging 
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sensors.  In  a  CCD  camera,  the  image  resolution  is  de¬ 
termined  by  the  size  of  its  photo-detector.  When  an  en¬ 
semble  of  several  shifted  images  are  available,  we  may 
reconstruct  an  image  with  higher  resolution  which  is 
equivalent  to  an  effective  increase  of  the  sampling  rate 
by  interpolation. 

The  problem  of  reconstructing  a  high  resolution  in¬ 
terpolated  image  from  a  sequence  of  undersarapled 
noisy  and  blurred  frames  has  been  tackled  by  proce¬ 
dures  developed  both  in  the  space  domain  and  the 
wavenumber  domain.  In  the  space  domain,  several 
methods  were  reported  [1,  2,  3].  Peleg,  ei  a/.,  [1]  pro¬ 
posed  an  interesting  scheme.  Starting  with  an  initial 
guess  for  a  high  resolution  image,  they  simulated  the 
imaging  process  to  compute  low  resolution  images.  The 
high  resolution  image  is  iteratively  improved  to  achieve 
simulated  low  resolution  images  closest  to  the  observed 
ones.  Tekalp,  et  al.^  [2]  proposed  a  scheme  based  on 
a  prior  POCS  (projection  onto  convex  sets)  method. 
Their  procedure  was  successfully  applied  to  blurred 
noisy  images.  However,  the  recursive  scheme  was  not 
explicit.  Ur  and  Gross  [3]  proposed  a  nonuniform  inter¬ 
polation  scheme  based  on  a  generalized  sampling  the¬ 
orem  to  obtain  an  improved  resolution  image  from  an 
ensemble  of  subpixel  level  shifted  low  resolution  images 
of  the  same  scene.  They  did  not  take  into  account  the 
presence  of  noise. 

Using  sequential  estimation  theory  in  the  wavenum¬ 
ber  domain,  an  efficient  method  was  developed  [4]  for 
recursively  updating  to  provide  satisfactory  reconstruc¬ 
tion  in  the  blur-free  case  provided  the  displacements  of 
the  frames  with  respect  to  a  reference  frame  were  either 
known  or  estimated.  It  was  observed  that  the  perfor¬ 
mance  deteriorated  when  the  blur  produced  zeros  in 
the  wavenumber  domain  and  theoretical  justification 
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for  this  can  be  provided.  The  problem  of  reconstruction 
in  the  wavenumber  domain  with  errors  present  both  in 
the  observation  and  data  was  tackled  by  the  method  of 
Total  Least  Squares  (TLS)  [5].  A  recursive  updating 
scheme  was  developed  which  performed  satisfactorily 
m  the  blur-free  case.  Very  recently,  Kim  [6]  proposed 
that  the  blurred  and  noisy  case  be  tackled  by  applying 
the  method  in  [5]  with  post-deblurring. 

In  this  paper,  we  expand  the  application  of  the  pro¬ 
cedures  in  [5]  and  [6]  to  a  wide  variety  of  blur  and  noise 
and  then  work  towards  a  theoretical  analysis  for  robust¬ 
ness  evaluation  of  the  TLS  method  in  the  wavenumber 
domain  where  the  data  is  complex-valued.  For  this  to 
be  possible,  the  results  from  the  TLS  theory  for  the 
case  of  real-valued  data  are  adapted  to  the  situation  at 
hand.  In  Section  2,  a  brief  review  of  high  resolution  im¬ 
age  reconstruction  by  the  TLS  algorithm  is  provided. 
Then,  the  robustness  of  the  TLS  algorithm  for  image 
reconstruction  is  evaluated  through  variance  analysis 
in  Section  3.  In  Section  4,  some  computer  simulation 
results  are  supplied  to  demonstrate  the  robustness  of 
the  TLS  algorithm  for  image  reconstruction.  Finally, 
concluding  remarks  are  made  in  Section  5. 

2.  HIGH  RESOLUTION  IMAGE 
RECONSTRUCTION  BY  THE  TLS 
ALGORITHM 

Suppose  that  the  Fourier  trcinsform,  F’*(«,  v)  of  the 
original  continuous  image  /(x,y)  is  appro4n’ated  by 
the  bandhmited  constraint 

IE  (u,  t;)|  =  0,  for  |u|  >  and  jvj  >  L^Wy,  (1) 

where  Lx,  Ly  are  some  finite  integers,  Wx  =  2itITx, 

Wy  ^J-KjTy,  and  T,,  Ty  are,  respectively,  the  sampling 
periods  along  the  a:  and  y  axes.  At  least  \LxLy  un- 
d^ampled  frames,  U{ix,iy),  i  =  1, 2, . .  .,4L,L„,  are 
needed  for  the  initial  reconstruction  of  f{x^y),  where 
and  iy  are  indices  of  natural  numbers.  To  accom¬ 
modate  explicitly  the  non-unity  sampling  periods  and 
shifts,  the  ith  frame  is  expressed  as 

fi{Xx,  ty)  =  f(ixTx  +  Sxi,iyTy  +  Syi) ,  (2) 

Where,,  =  0,1,. =  0,1 . Y  -  1,  and 

<)«.•,  6yi  are,  respectively,  the  estimated  shifts  along  the 
a;  and  y  axes.  It  is  assumed  that  diflFerences  between 
the  true  shifts  and  the  estimated  ones  along  the  x  and 
y  «es  are,  respectively,  A4.  and  AJ,,-.  If  there  are 
k  frames  a,vailable,  the  multiframe  image  restoration  , 
model  IS  [5] 

Z=[^  +  E]F-i-N,  (3)  i 


where 


[Zi  Z2  ...  ZkY,  #  =  [Yi  Yj 

[4>il  <{>i2  ■  .  . 

1  /  /  /  m  / 


i  (s 


(r-l)mod(2i,)-L^, 
L(r-l)/(2L,)J-£^, 
[s’!!  Si2  .  .  ■  SipY, 

iO.TT  ( \X  .  ( 


=  j2ir  (  AJ, 


lx 

MTx  Tx 


E  =  [AYi  AY2  ...  AY*]', 

N  =  [N1N2  ...  Nkf, 

^  =  [^mf.(l)  E7.„(2)  ...  Fm„(p)]‘. 

Furthermore^  the  superscript  t  denotes  the  transpose 

operator,  p  _  ALxLy,  Zi  and  AT,-  are,  respectively,  the 
values  of  the  DFTs  of  the  ,‘th  noisy  undersampled  frame 
Md  the  additive  noise  at  a  generic  index  2-tuple  (m  n] 
in  the  wavenumber  domain.  is  the  ith  interpcv 

lated  component  at  index  2-tuple  (m,  n).  This  2-tuple 
IS  suppressed  in  Z,-,  Ni,  Z  and  N  for  brevity.  Without 
loss  of  generality,  A,-  is  assumed  to  be  zero  mean. 

In  [5,  6],  the  TLS  algorithm  [7]  developed  for  the 
real-valued  data  is  successfully  adapted  for  high  reso¬ 
lution  image  reconstruction  using  the  model  given  in 
Equation  (3). 


3.  PERFORMANCE  ANALYSIS  OF  THE 
TLS  ALGORITHM  FOR  IMAGE 
RECONSTRUCTION 

The  TLS  theory  has  been  studied  extensively  in  the 
ca^  of  real-valued  data  [7].  In  [5]  and  [6],  the  complex- 
va  ued  data  problem  has  been  transformed  to  an  equiv¬ 
alent  real-valued  data  problem  to  which  the  TLS  the¬ 
ory  is  applied.  In  this  paper,  theoretical  analysis  of  the 
performance  of  the  TLS  algorithm  [5,  6]  for  complex- 
valued  data  is  proposed.  It  is  shown  that  the  errors 
originating  from  the  additive  noise  Ni  and  the  inaccu¬ 
rate  estimates  oiSxi  and  Syi  can  be  suppressed  simulta¬ 
neously  by  applying  the  TLS  algorithm  [5,  6],  provided 
certain  constraints  are  satisfied. 

Let  Cr  and  Ci  denote,  respectively,  the  real  and 
imaginary  parts  of  a  complex  number  C.  Then,  (3) 
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I  be  rewritten  as  (4)  and  (5)  below. 

2,  =  [^r+Er  -^i  -  +  NrjS) 

Zi  =  [^i  +  Ei  ^  +  Ni.  (5) 

pfbe  preceding  two  equations  may  be  combined  as 

U  =  [A  +  AA]V  +  A»r.  (6) 


AX  1 

Ei  Er  )  ■ 


i  Suppose  that  there  are  k  noisy  undersampled  images 
■  available.  It  has  been  shown  [7,  p.  242]  that 

|:-  coy{V)!=^{l  +  \\Vf)<T^A*A-ka^l)-\  (7) 

|if  the  row  vectors  of  [AA  AW]  are  independent  and 
|identic2illy  distributed  with  common  zero  mean  vector 
flkid  common  covariance  matrix  <t^I  where  (>  0)  is 
^a  small  unknown  constant.  In  (7), 

M- 

)'  V 

|Equation  (8)  can  be  further  simplified  by  making  use 
|pf  (9)  and  (10)  below.  As  the  number  of  frames  k  in- 
l^creases,  it  is  not  difficult  to  verify  that  the  Hermitian 
matrix  where  the  superscript  H  denotes  Herini- 

tian  conjugate,  asymptotically  approaches  a  diagonal 
-matrix  with  nonnegative  entries.  Therefore,  for  a  large 
ilbut  finite  fc, 


\^x^y 

'On  the  other  hand, 

-  ^i^r) 

From(9)  and  (10),  we  C3m  have 


(9) 


(10) 


r  + 


diag 


•••  TITJ, 


^pxp 

(11a) 


0. 


(11b) 


Then,  from  (8),  (Ha),  and  (11b),  A‘A  approximates 
a  2p  X  2p  diagonal  matrix  as  below. 


a‘a 


diag 


7^  7^2 


rr2rr2 

-^x^y 


2px2p 


(12) 


By  making  use  of  (12),  (7)  can  be  rewritten  as 
1  cr^T^T^ 

cov(V)=,j(l  +  imP)5-^X.  (U) 

If  the  additive  noise  is  Gaussian  and  is  small,  then 
it  has  been  shown  that  cov(V)  can  reach  the  Cramer- 
Rao  lower  bound  [7,  p.  243].  In  other  words,  the  image 
reconstructed  using  the  TLS  algorithm  [5,  6]  has  min¬ 
imum  variance  with  respect  to  all  unbiased  estimates. 
As  the  number  of  undersampled  images,  increases, 
the  quality  of  the  reconstructed  image  becomes  bet¬ 
ter  and  better  as  shown  in  (13).  In  the  next  section, 
computer  simulation  results  are  provided. 


4.  COMPUTER  SIMULATION  RESULTS 

In  this  computer-simulated  example,  a  high  resolution 
256  X  256  image,  as  shown  in  Figure  1,  is  recursively 
reconstructed  from  a  set  of  low  resolution  128  x  128 
noisy  frames  which  are  shifted  with  respect  to  a  ref¬ 
erence  frame.  These  shifts  or  displacements  are  not 
accurately  known.  In  order  to  reduce  the  size  of  the 
arrays  to  be  processed,  the  input  image  was  partitioned 
into  16  nonoverlapping  sections  each  of  size  32  x  32  and 
the  recursive  TLS  algorithm  developed  in  [5,  6]  is  inde¬ 
pendently  applied  to  each  such  section.  For  each  one 
of  these  16  sections,  the  interpolation  problem  corre¬ 
sponds  to  the  reconstruction  of  a  64  x  64  image  from 
a  sequence  of  shifted  low  resolution  32  x  32  noisy  in¬ 
put  frames  when  the  interframe  displacements  are  not 
accurately  known.  To  generate  the  k  =  16  shifted  low 
resolution  input  frames  required  in  the  simulation  from 
the  available  data,  we  use  the  DFT-based  interpolation 
technique  described  in  [6].  After  assigning  one  of  the 
input  frames  to  be  the  reference  frame,  we  label  it  as 
frame  number  1,  and  the  remaining  ones  are  sequen¬ 
tially  labeled  from  i  =  2  to  t  =  16.  The  relative  shifts  of 
these  frames  with  respect  to  the  frame  number  1  along 
the  X  and  y  axes  are  denoted,  respectively,  by  S^i  and 
Syi  for  i  =  2, 3, . . . ,  16.  The  estimation  errors,  AJ**  and 
ASyi,  of  5^i  and  Syi,  i  =  2, 3, . . . ,  16,  are  assumed  to  be 
uniformly  distributed  over  [-|,  |]*  Subsequently,  each 
frame  is  corrupted  by  additive  noise  with  SNR  level  of 
20  dB.  To  illustrate  the  intermediate  steps  of  the  recon¬ 
struction,  we  provide  a  generic  set  of  sixteen  32  x  32 
frames,  as  shown  in  Figure  2,  corresponding  to  a  sec¬ 
tion  of  the  low  resolution  noisy  image.  From  this  set,  a 
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64  X  64  image  is  obtained  by  application  of  the  recur¬ 
sive  TLS  algorithm.  The  sequence  of  estimates  leading 
to  the  construction  of  the  high  resolution  filtered  im¬ 
age  is  shown  in  Figure  3.  As  expected,  the  first  few 
iterations  of  the  recursive  TLS  algorithm  provide  poor 
estimates  because  the  algorithm  is  in  a  transient  stage. 
Subsequently,  the  estimates  improve. 

When  the  undersampled  images  are  blurred,  a  post¬ 
deblurring  approach  using  Wiener  filtering  is  proposed 
in  [6].  In  this  computer  simulation,  the  frames  are  cor¬ 
rupted  by  a  6-pixel  uniform  linear  motion  blur  and  30 
dB  SNR  additive  noise.  Figures  4  and  5  show,  respec¬ 
tively,  the  high  resolution  images  obtained  by  applying 
the  recursive  TLS  algorithm  without  post-deblurring 
and  with  post-deblurring. 

5.  CONCLUDING  REMARKS 

In  computer  simulations,  the  resolution  of  the  recon¬ 
structed  image  increased  noticeably  when  the  number 
of  noisy  undersampled  frames  reaches  ALj:Ly  (=:  p),  the 
minimum  number  of  undersampled  frames  needed.  As 
k  increases,  the  quality  of  the  reconstructed  image  con¬ 
tinues  to  improve.  However,  the  improvement  reaches 
a  level  of  saturation  after  several  frames  k  >  p.  This 
may  be  due  to  the  assumption  made  to  derive  (7). 

As  pointed  out  in  [6],  the  TLS  algorithm  outper¬ 
forms  the  constrained  least  squares  algorithm  when  the 
blurs  in  the  undersampled  frames  au-e  identical.  If  the 
undersampled  frames  are  obtained  from  the  same  cam¬ 
era,  it  is  reasonable  to  assume  that  the  blurs  in  dif¬ 
ferent  frames  are  identical.  The  performance  analy¬ 
sis  proposed  here  further  demonstrates  the  advantage 
of  using  the  TLS  algorithm  in  high  resolution  image 
reconstruction  from  noisy  undersampled  frames  when 
the  displacement  of  any  frame  with  respect  to  a  ref¬ 
erence  frame  is  not  accurately  estimated.  Extensive 
simulations  conducted  and  reported  in  [6]  support  the 
implications  of  performance  analysis  reported  in  this 
paper.  Future  research  is  planned  on  generalizing  the 
techniques  in  [5],  [6],  and  this  paper  to  the  case  of  mul- 
tispectral  frames. 
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•e  2:  Sixteen  32  x  32  input  girl  frames 


Figure  4;  A  64  x  64  girl  image  reconstruction  sequence 
without  post-deblurring 


Figure  3;  A  64  x  64  girl  image  reconstruction  sequence 


Figure  5:  A  64  x  64  girl  image  reconstruction  sequence 
with  post-deblurring 


